{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling the potential energy surface "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction \n",
    "\n",
    "This interactive notebook demonstrates how to utilize the Potential Energy Surface (PES) samplers algorithm of Qiskit Nature to generate the dissociation profile of a molecule. We use the Born-Oppenheimer Potential Energy Surface (BOPES) and demonstrate how to exploit bootstrapping and extrapolation to reduce the total number of function evaluations in computing the PES using the Variational Quantum Eigensolver (VQE). \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import common packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "\n",
    "# qiskit\n",
    "from qiskit.utils import QuantumInstance, algorithm_globals\n",
    "from qiskit import Aer\n",
    "from qiskit.algorithms.minimum_eigen_solvers import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.algorithms.optimizers import SLSQP\n",
    "from qiskit.circuit.library import ExcitationPreserving\n",
    "from qiskit import BasicAer\n",
    "from qiskit.algorithms import NumPyMinimumEigensolver, VQE\n",
    "from qiskit.algorithms.optimizers import SLSQP\n",
    "\n",
    "# qiskit nature imports\n",
    "from qiskit_nature.problems.second_quantization import ElectronicStructureProblem\n",
    "from qiskit_nature.converters.second_quantization import QubitConverter\n",
    "from qiskit_nature.mappers.second_quantization import JordanWignerMapper\n",
    "from qiskit_nature.algorithms import GroundStateEigensolver\n",
    "from qiskit_nature.drivers import UnitsType, Molecule\n",
    "from qiskit_nature.drivers.second_quantization import (\n",
    "    ElectronicStructureDriverType,\n",
    "    ElectronicStructureMoleculeDriver,\n",
    ")\n",
    "from qiskit_nature.algorithms.pes_samplers import BOPESSampler, Extrapolator\n",
    "\n",
    "import warnings\n",
    "\n",
    "warnings.simplefilter(\"ignore\", np.RankWarning)\n",
    "\n",
    "algorithm_globals.random_seed = 75"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we use the H2 molecule as a model system for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "stretch1 = partial(Molecule.absolute_stretching, atom_pair=(1, 0))\n",
    "mol = Molecule(\n",
    "    geometry=[(\"H\", [0.0, 0.0, 0.0]), (\"H\", [0.0, 0.0, 0.3])],\n",
    "    degrees_of_freedom=[stretch1],\n",
    ")\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = ElectronicStructureMoleculeDriver(mol, driver_type=ElectronicStructureDriverType.PYSCF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.3])]\n"
     ]
    }
   ],
   "source": [
    "print(mol.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a perturbation to the molecule along the absolute_stretching dof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.5])]\n",
      "[('H', [0.0, 0.0, 0.0]), ('H', [0.0, 0.0, 0.8999999999999999])]\n"
     ]
    }
   ],
   "source": [
    "mol.perturbations = [0.2]\n",
    "print(mol.geometry)\n",
    "\n",
    "mol.perturbations = [0.6]\n",
    "print(mol.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate bond dissociation profile using BOPES Sampler\n",
    "\n",
    "Here, we pass the molecular information and the VQE to a built-in type called the BOPES Sampler. The BOPES Sampler allows the computation of the potential energy surface for a specified set of degrees of freedom/points of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First we compare no bootstrapping vs bootstrapping \n",
    "\n",
    "Bootstrapping the BOPES sampler involves utilizing the optimal variational parameters for a given degree of freedom, say r (ex. interatomic distance) as the initial point for VQE at a later degree of freedom, say r + $\\epsilon$. By default, if boostrapping is set to True, all previous optimal parameters are used as initial points for the next runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "es_problem = ElectronicStructureProblem(driver)\n",
    "\n",
    "qubit_converter = QubitConverter(JordanWignerMapper())\n",
    "quantum_instance = QuantumInstance(backend=Aer.get_backend(\"aer_simulator_statevector\"))\n",
    "solver = VQE(quantum_instance=quantum_instance)\n",
    "\n",
    "me_gsc = GroundStateEigensolver(qubit_converter, solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance1 = partial(Molecule.absolute_distance, atom_pair=(1, 0))\n",
    "mol = Molecule(\n",
    "    geometry=[(\"H\", [0.0, 0.0, 0.0]), (\"H\", [0.0, 0.0, 0.3])],\n",
    "    degrees_of_freedom=[distance1],\n",
    ")\n",
    "\n",
    "# pass molecule to PSYCF driver\n",
    "driver = ElectronicStructureMoleculeDriver(mol, driver_type=ElectronicStructureDriverType.PYSCF)\n",
    "\n",
    "es_problem = ElectronicStructureProblem(driver)\n",
    "\n",
    "\n",
    "# Specify degree of freedom (points of interest)\n",
    "points = np.linspace(0.25, 2, 30)\n",
    "results_full = {}  # full dictionary of results for each condition\n",
    "results = {}  # dictionary of (point,energy) results for each condition\n",
    "conditions = {False: \"no bootstrapping\", True: \"bootstrapping\"}\n",
    "\n",
    "\n",
    "for value, bootstrap in conditions.items():\n",
    "    # define instance to sampler\n",
    "    bs = BOPESSampler(gss=me_gsc, bootstrap=value, num_bootstrap=None, extrapolator=None)\n",
    "    # execute\n",
    "    res = bs.sample(es_problem, points)\n",
    "    results_full[f\"{bootstrap}\"] = res.raw_results\n",
    "    results[f\"points_{bootstrap}\"] = res.points\n",
    "    results[f\"energies_{bootstrap}\"] = res.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare to classical eigensolver\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define numpy solver\n",
    "solver_numpy = NumPyMinimumEigensolver()\n",
    "me_gsc_numpy = GroundStateEigensolver(qubit_converter, solver_numpy)\n",
    "bs_classical = BOPESSampler(\n",
    "    gss=me_gsc_numpy, bootstrap=False, num_bootstrap=None, extrapolator=None\n",
    ")\n",
    "# execute\n",
    "res_np = bs_classical.sample(es_problem, points)\n",
    "results_full[\"np\"] = res_np.raw_results\n",
    "results[\"points_np\"] = res_np.points\n",
    "results[\"energies_np\"] = res_np.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Energy')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4VOX5//H3PZONLGRnCUsW9iUhQNhVFtG6KwJWKxWkLmir1X5t1dqvX639tdhate51Rawibqi1KBUEEUEg7GHfAoQ9CQkJIcvMPL8/5hBDSMgAmZwk3K/rmosz55w55zMnw9xztucRYwxKKaWULxx2B1BKKdV0aNFQSinlMy0aSimlfKZFQymllM+0aCillPKZFg2llFI+06KhmhQReUVE/teG9X4pIhPP8rUdRaRYRJz1nau+icgYEdlj5e0rIutFZIQ17TER+ZfNEZXNRO/TUI2FiGQDrQEX4AY2ANOBV40xHhujnRHrfdxmjJlrd5YzJSLbgd8YYz6rYdpjQGdjzIQGD6YaDd3TUI3N1caYCCARmAo8CLxhb6TmQUQCfJgtEVjv7yyq6dKioRolY0yhMeZz4KfARBHpDSAi00TkT9ZwnIh8ISIFIpIvIt+JiMOa9qCI7BWRIhHZLCIXW+ODReRZEdlnPZ4VkeAT6xWRa0VktYgcFZHtInKZNX6BiNxmDXcSkW9EJE9EckXkXRGJsqa9A3QE/m0d4vmdiCSJiDnxpS0iCSLyuZV5m4jcXmX9j4nIByIy3cq+XkQyattO1nLvFZEdVpa/VdkGk0TkexF5RkTygMdExCEifxCRXSJyyFpPpLVdigEnsMba40BEskVkdC3rHiwii63tv+bEYSzVvGnRUI2aMWYZkANcWMPk/7GmxeM9rPV7wIhIN+BXwABrr+UnQLb1mkeAwUA60AcYCPwBQEQG4j0c9lsgCrioyuuqEuAvQALQA+gAPGbl/TmwG+8eU7gx5q81vP59K3cCMA74s4iMqjL9GmueKOBz4IWat06lMUAG0A+4FphcZdogYAfe7fP/gEnWYySQAoQDLxhjyowx4dZr+hhjOp1uhSLSDvgP8CcgBngA+FhE4uvIqpo4LRqqKdiH94upugqgLZBojKkwxnxnvCfp3EAw0FNEAo0x2caY7dZrbgb+aIw5ZIw5DDwO/Nya9gvgTWPM18YYjzFmrzFmU/WVGmO2WfOUWct4GhjuyxsRkQ7AMOBBY0ypMWY18DpwS5XZFhljZhtj3MA7eIvb6TxpjMk3xuwGngVuqjJtnzHmeWOMyxhz3Hr/TxtjdhhjioGHgRt9PHRV1QRgtpXTY4z5GsgErjjD5agmRouGagraAfk1jP8bsA34r3V45iHwfqkD9+H99X9IRN4XkQTrNQnArirL2GWNA+8ew3bqICKtrWXuFZGjwL+AOB/fSwKQb4wpqpahXZXnB6oMlwAhdXyp76m2rIRapp1Yf/X3H4B3T+RMJALjrUNTBSJSAFyAt4irZkyLhmrURGQA3i/URdWnGWOKjDH/Y4xJwXtI5zcnzl0YY94zxlyA98vNAE9aL9tnjTuhozUOvF+wpz0sY/mztcxUY0xLvL+6pWq007x2HxAjIhHVMuz1Yb216VBtWfuqPK+epab37wIOnuE69wDvGGOiqjzCjDFTz3A5qonRoqEaJRFpKSJX4T22/y9jzLoa5rlKRDqLiACFeA9LeUSkm4iMsk5wlwLHgROX7M4A/iAi8SISBzyKd08BvFdp3SoiF1snjNuJSPca4kUAxUChdWz/t9WmH8R7vuAUxpg9wGLgLyISIiJpeA+Lncv9D78VkWjr0NevgZmnmXcGcL+IJItION4CONMY4zrDdf4LuFpEfiIiTuu9jBCR9mf3FlRToUVDNTb/FpEivL9kH8F7vuDWWubtAszF+wW+BHjJGDMf7/mMqUAu3kM9rfAeuwfvidtMYC2wDlhpjTtx0v1W4Bm8RehbTv5VfsLjeE86F+I9GfxJtel/wVuYCkTkgRpefxOQhPdX/yzg/87xno7PgBXAaivP6S5RfhPveZKFwE68RfWeM12hVfyuxXvxwWG8f6/fot8pzZ7e3KdUEyYiBuhincdRyu/0V4FSSimfadFQSinlMz08pZRSyme6p6GUUspnZ3oXaKMXFxdnkpKS7I6hlFJNyooVK3KNMXU2A9PsikZSUhKZmZl2x1BKqSZFRHbVPZcenlJKKXUGtGgopZTymRYNpZRSPmt25zSUUuemoqKCnJwcSktL7Y6i/CAkJIT27dsTGBh4Vq/XoqGUOklOTg4REREkJSXhbQtSNRfGGPLy8sjJySE5OfmslqGHp5RSJyktLSU2NlYLRjMkIsTGxp7TXqQWDaXUKbRgNF/n+rfVomHJOZTNI2+NYfai6XZHUUqpRkuLhsUhwueObSzbPtvuKEqpczBp0iQ++uijc17O6tWrmT379N8H2dnZvPfee+e8Ll+88sorTJ9u/49aLRqWhPhEWro95Jafaa+XSqnm6FyLhst1pp0hnt6UKVO45ZZb6nWZZ0OLRhXxLidHPEftjqHUeS07O5sePXpw++2306tXLy699FKOHz8OeL/IBw8eTFpaGmPGjOHIkSM1LmPu3LlkZGTQtWtXvvjiC8B7gv/WW28lNTWVvn37Mn/+/FrHl5eX8+ijjzJz5kzS09OZOXMm3377Lenp6aSnp9O3b1+Kiop46KGH+O6770hPT+eZZ55h2rRpXHPNNYwaNYqLL76Y4uJiLr74Yvr160dqaiqfffZZ5Xvs3r07N998Mz169GDcuHGUlJQA3qaQfve735GamsrAgQPZts3bv9Zjjz3GU089BcCIESN48MEHGThwIF27duW7774DoKSkhBtuuIGePXsyZswYBg0aVO/NKuklt1XEmFD2OovtjqFUo/H4v9ezYV/9/pDqmdCS/7u612nn2bp1KzNmzOC1117jhhtu4OOPP2bChAnccsstPP/88wwfPpxHH32Uxx9/nGefffaU12dnZ7Ns2TK2b9/OyJEj2bZtGy+++CIiwrp169i0aROXXnopW7ZsqXX8H//4RzIzM3nhhRcAuPrqq3nxxRcZNmwYxcXFhISEMHXqVJ566qnKwjRt2jRWrlzJ2rVriYmJweVyMWvWLFq2bElubi6DBw/mmmuuAWDz5s288cYbDBs2jMmTJ/PSSy/xwAPe3oEjIyNZt24d06dP57777qtcflUul4tly5Yxe/ZsHn/8cebOnctLL71EdHQ0GzZsICsri/T09HP6W9VE9zSqiAmI5nAAlJeX2R1FqfNacnJy5Rde//79yc7OprCwkIKCAoYPHw7AxIkTWbhwYY2vv+GGG3A4HHTp0oWUlBQ2bdrEokWLmDBhAgDdu3cnMTGRLVu21Dq+umHDhvGb3/yG5557joKCAgICav7NfckllxATEwN474v4/e9/T1paGqNHj2bv3r0cPOg9BN6hQweGDRsGwIQJE1i0aFHlMm666abKf5csWVLjeq6//vqTtg/AokWLuPHGGwHo3bs3aWlpNb72XOieRhWtWrSnoiKHTbtWkdZlsN1xlLJdXXsE/hIcHFw57HQ6Kw9P+ar6ZaX1cQnxQw89xJVXXsns2bMZNmwYc+bMqXG+sLCwyuF3332Xw4cPs2LFCgIDA0lKSqq8R+J0GWsbrurENnI6nfV+/uR0dE+jinYxXQDYsnuFzUmUUtVFRkYSHR1defz+nXfeqdzrqO7DDz/E4/Gwfft2duzYQbdu3bjwwgt59913AdiyZQu7d+8+7fiIiAiKiooql7l9+3ZSU1N58MEHGTBgAJs2bTplnuoKCwtp1aoVgYGBzJ8/n127fmx9fPfu3ZV7Ee+99x4XXHBB5bSZM2dW/jtkyBCft9GwYcP44IMPANiwYQPr1q3z+bW+0j2NKrq06wcH32F37ka7oyilavD2228zZcoUSkpKSElJ4a233qpxvo4dOzJw4ECOHj3KK6+8QkhICHfffTd33XUXqampBAQEMG3aNIKDg2sdP3LkSKZOnUp6ejoPP/wwixYtYv78+TgcDnr16sXll1+Ow+HA6XTSp08fJk2aRHR09Ek5br75Zq6++mpSU1PJyMige/fuldO6devGiy++yOTJk+nZsyd33XVX5bQjR46QlpZGcHAwM2bM8Hn73H333UycOJGePXvSvXt3evXqRWRk5Blu5dNrdn2EZ2RkmLO9WuBYSRFDPhjClaYTf7n1s3pOplTTsHHjRnr06GF3jGYtOzubq666iqysrFOmnehILi4u7oyX63a7qaioICQkhO3btzN69Gg2b95MUFDQSfPV9DcWkRXGmIy61qF7GlWEhUYQ7zbkuXPtjqKUUmespKSEkSNHUlFRgTGGl1566ZSCca60aFQT6w4kH73sVinlP0lJSTXuZQCVV0KdjYiICL93d60nwquJIYJcZ8NdiaCUUk2JFo1qYgPjyQtwUFCkh6iUUqo6LRrVtAlPBCBr2w82J1FKqcZHi0Y1HeO9VxRs37fG5iRKKdX4aNGopkfyQAD2FWyzOYlS56/s7Gx69+59zst59tlnKxsCrM20adPYt2/fOa/LF0OHDm2Q9fiTFo1qOrdPJdhjOHy8YT5ESin/Odei4Xa76zXP4sWL63V5dtCiUY3D6aSVC/LdNTe5rJRqGC6X65Smw+fNm0ffvn1JTU1l8uTJlJV5Gxetafxzzz3Hvn37GDlyJCNHjsTtdjNp0iR69+5NamoqzzzzDB999BGZmZncfPPNpKenc/z4cZKSknjwwQfp168fH374Ia+99hoDBgygT58+jB07trIITZo0iSlTppzSBPu0adO49tprGTFiBF26dOHxxx+vfE/h4eEALFiwgBEjRjBu3LjKJtJP3Gg9e/ZsunfvTv/+/bn33nu56qqrGnKz18mW+zREJAaYCSQB2cANxpgav6VFpCWwAfjUGPOrhsgX6wkhz3FmDaQp1Sx9+RAcqOf2i9qkwuVT65ytetPhTz/9NP/85z+ZN28eXbt25ZZbbuHll19mypQpTJo06ZTx9913H08//TTz588nLi6OFStWsHfv3sr7IwoKCoiKiuKFF17gqaeeIiPjx5uhY2NjWblyJQB5eXncfvvtAPzhD3/gjTfe4J577gFqboIdYNmyZWRlZREaGsqAAQO48sorT1o+wKpVq1i/fj0JCQkMGzaM77//noyMDO68804WLlxIcnJyZWu3jYldexoPAfOMMV2Aedbz2jwB1Nz+sZ/EOCI5HODBU8+7pkop31VvOnzevHkkJyfTtWtX4Mem0Tdv3lzj+OpSUlLYsWMH99xzD1999RUtW7asdd0//elPK4ezsrK48MILSU1N5d1332X9+vWV02pqgh28zaPHxsbSokULrr/++pOaPT9h4MCBtG/fHofDQXp6OtnZ2WzatImUlBSSk5MBGmXRsOuO8GuBEdbw28AC4MHqM4lIf6A18BVQZ5so9SUupC3HPIfYc3A7iQldG2q1SjU+PuwR+Ev1JsGjoqLIy8s76+VFR0ezZs0a5syZwyuvvMIHH3zAm2++WeO8VZs3nzRpEp9++il9+vRh2rRpLFiwoNaMJ5770jR79ebfG7J583Nh155Ga2PMfmv4AN7CcBIRcQB/Bx6oa2EicoeIZIpI5uHDh885XNvIFAA27Fx6zstSSp2d6k2HZ2RkkJ2dXXkI6ETT6N26datxPHBS0+W5ubl4PB7Gjh3Ln/70p8rDT3U1b15UVETbtm2pqKiobEL9hJqaYAf4+uuvyc/P5/jx43z66aeVe0x16datGzt27KhsSuREE+mNid/2NERkLtCmhkmPVH1ijDEiUlNTu3cDs40xOXV1oGKMeRV4Fbyt3J5d4h+ltEmDI7PIPlhz2zBKKf+r3nT4c889x+DBgxk/fjwul4sBAwYwZcoUgoODeeutt04ZD3DHHXdw2WWXkZCQwLPPPsutt96Kx+MB4C9/+Qvw4wntFi1a1NhL3hNPPMGgQYOIj49n0KBBJxWYmppgB++hp7Fjx5KTk8OECRNOOZ9RmxYtWvDSSy9x2WWXERYWxoABA85pG/qDLU2ji8hmYIQxZr+ItAUWGGO6VZvnXeBCwAOEA0HAS8aY053/OKem0U84mLeX0V9cxvX05PGJja/SK+VP2jS6byZNmsRVV13FuHHjTho/bdq0k/oWP1PFxcWEh4djjOGXv/wlXbp04f7776+PyJXOpWl0uw5PfQ5MtIYnAqd0XmGMudkY09EYk4T3ENX0ugpGfWkd245It4fcsoMNsTqllKr02muvkZ6eTq9evSgsLOTOO++0O9JJ7NrTiAU+ADoCu/BecpsvIhnAFGPMbdXmnwRk+HLJbX3saQCMeTWNFiaA9+5cec7LUqop0T2N5q/JdcJkjMkDLq5hfCZwWw3jpwHT/B6sihjC2OOs/eSYUkqdj/SO8FrEOGM4HADl5WV2R1FKqUZDi0YtWoW2xyXChp3L7Y6ilFKNhhaNWrSL8d7Ut2W3ntNQSqkTtGjUomv7/gDk5G+2OYlSSjUeWjRq0bPTABzGcOjYHrujKKVUo6FFoxahIWHEuwx5Lu0rXKmGlp2dTY8ePbj99tvp1asXl156KcePH2fEiBGcuKQ+NzeXpKQkwHtD3XXXXccll1xCUlISL7zwAk8//TR9+/Zl8ODB5OfnAzBixAh+/etfk56eTu/evVm2bBkej4cuXbpwogkij8dD586dqY8miZojuxosbBLiPEHkc8zuGErZ5sllT7Ipf1O9LrN7THceHHhK+6Sn2Lp1KzNmzOC1117jhhtu4OOPPz7t/FlZWaxatYrS0lI6d+7Mk08+yapVq7j//vuZPn069913HwAlJSWsXr2ahQsXMnnyZLKyspgwYQLvvvsu9913H3PnzqVPnz7Ex8fXy/ttbnRP4zRiiCDX2TRanlSquUlOTiY9PR2A/v37VzbiV5uRI0cSERFBfHw8kZGRXH311QCkpqae9NoTzY1fdNFFHD16lIKCAiZPnsz06dMBePPNN7n11lvr/w01E7qncRoxQfHkyxGOFB4mOlJ/dajzjy97BP5Svenw48ePExAQUNngYGlpaa3zOxyOyucOh+OkZsdrara8Q4cOtG7dmm+++YZly5ad0pqt+pHuaZxGmwhvRyhZ209t+VIp1fCSkpJYsWIFAB999NFZLeNEc+OLFi0iMjKSyMhIAG677TYmTJjA+PHjcTqd9RO4GdKicRod471ts2zbt8bmJEopgAceeICXX36Zvn37kpt7dhephISE0LdvX6ZMmcIbb7xROf6aa66huLhYD03VwZYGC/2pvhosBNi2O4sx82/iRmd/HpkwrV6WqVRj15wbLBwxYsQp/YGfkJmZyf333893331nQ7KG1eQaLGwqUtr1IMRjyK3YZ3cUpZQfTZ06lZdfflnPZfhAi8ZpOJxOWrmEfM8Ru6MopepB1f69q3rooYd46KEG6a6nydNzGnWINSHkOUrrnlGpZqS5HbZWPzrXv60WjTrESCSHAzx43G67oyjVIEJCQsjLy9PC0QwZY8jLy6vsy/xs6OGpOsSFJFDiOciuA1tJbtfd7jhK+V379u3JycnRZjSaqZCQENq3b3/Wr9eiUYeEqE6Qv4oNO3/QoqHOC4GBgSQnJ9sdQzVSeniqDilt0wDYdWiDzUmUUsp+WjTq0KvTIAAOFO20OYlSStlPD0/VIT46gWi3hzz3IbujKKWU7bRo+CDO5SSfo3bHUEop2+nhKR/EmDDyHBV2x1BKKdtp0fBBTGAshwKgtKzE7ihKKWUrLRo+aB3aAbcIG3ausDuKUkrZSouGD9rFdAVgy+76aT1XKaWaKi0aPuja0dta8N78LTYnUUope2nR8EHP5P44jeFgyR67oyillK30klsfhASH0soF+e48u6MopZStdE/DR7GeQPLlmN0xlFLKVlo0fBQjLckN0ObRlVLnNy0aPooNasURp4O8ggN2R1FKKdto0fBRmwhvU9Hrti2xOYlSStlHi4aPElv1BGDH/nU2J1FKKfto0fBRj8QBAOwr2GZzEqWUso9ecuujpITuhHo85FbstzuKUkrZRouGjxxOJ/EuB/mmwO4oSillGz08dQZiPSHkS6ndMZRSyjZaNM5AjDOagwEGj1vv11BKnZ9sKRoiEiMiX4vIVuvf6Frm6ygi/xWRjSKyQUSSGjbpyeJCEih1CDv2brQzhlJK2cauPY2HgHnGmC7APOt5TaYDfzPG9AAGArZ21J0Q1QmADdlL7YyhlFK2satoXAu8bQ2/DVxXfQYR6QkEGGO+BjDGFBtjbO06r3NCHwB2H9Y9DaXU+cmuotHaGHPi2tUDQOsa5ukKFIjIJyKySkT+JiLOmhYmIneISKaIZB4+fNhfmendaYg3cNFOv61DKaUaM79dcisic4E2NUx6pOoTY4wREVPDfAHAhUBfYDcwE5gEvFF9RmPMq8CrABkZGTUtq15ER8YT4/KQ7/ZfYVJKqcbMb0XDGDO6tmkiclBE2hpj9otIW2o+V5EDrDbG7LBe8ykwmBqKRkOKcweQT5GdEZRSyjZ2HZ76HJhoDU8EPqthnuVAlIjEW89HARsaINtpxRBGrqPc7hhKKWULu4rGVOASEdkKjLaeIyIZIvI6gDHGDTwAzBORdYAAr9mUt1JsQBy5AUJJqXbIpJQ6/9jSjIgxJg+4uIbxmcBtVZ5/DaQ1YLQ6tQrrgLt0Jxu2Lyej1wi74yilVIPSO8LPUPuYbgBs27vS5iRKKdXwtGicoe6JGQDsydtscxKllGp4WjTOUPek/gQYw+GSHLujKKVUg9Om0c9QUFAwrVyQ58m3O4pSSjU43dM4C7HuII6gV08ppc4/WjTOQowjkgMBbm0iXSl13tGicRaSI7pT5HSwbMM8u6MopVSD0qJxFjK6XAbA0k3/sTmJUko1LC0aZ2FYnysI83jYVrDW7ihKKdWgfCoaVvPkV4qIFhkgICCQ5IpgdpNrdxSllGpQvhaBl4CfAVtFZKqIdPNjpiYhMaA9uwINeQUH7I6ilFINxqeiYYyZa4y5GegHZANzRWSxiNwqIoH+DNhYdW89ALcI81d8bHcUpZRqMD4fbhKRWLydIN0GrAL+gbeIfO2XZI3c8PTxAGTlfGdzEqWUajg+3REuIrOAbsA7wNVVumqdKSKZ/grXmCW3605ChSHbvcPuKEop1WB8bUbkOWPM/JomGGMy6jFPk5LoiWRzQAEetxuHs8buy5VSqlnxtWhEi8j11cYVAuuMMTV11XpeSAnrxhLXctZuW0J6twvsjqOUUn7n6zmNXwCvAzdbj9eAB4HvReTnfsrW6PVN8XaDvnj95zYnUUqphuFr0QgEehhjxhpjxgI9AQMMwls8zksX9r2WYI9ha/5qu6MopVSD8LVotDfGHKzy/BDQwRiTD1TUf6ymITQkjOSKAHZ7DtY9s1JKNQO+ntNYICJfAB9az8da48KAAr8kayISHW35JmAPRccKiAiLsjuOUkr5la97Gr8E3gLSrcd04JfGmGPGmJH+CtcUdI3PoEJv8lNKnSfqLBoi4gS+McZ8bIy533p8ZIwxDZCv0bsw9ToA1u5aYG8QpZRqAHUWDWOMG/CISGQD5GlyeqT0p5XLw87j2+yOopRSfufrOY1iYJ2IfA0/9nNqjLnXL6mamCRXBLucR+2OoZRSfudr0fjEeqgaJIV2ZplnDZt3rqJbcl+74yillN/4VDSMMW+LSAugozFms58zNTlpHUfwQfYavls3S4uGUqpZ87UTpquB1cBX1vN0EdHboC0j+l9PgDFsOrzc7ihKKeVXvl5y+xgwEOueDGPMaiDFT5manMjwGJLLHex27697ZqWUasJ8LRoVxpjCauM89R2mKevoaM2OQBelZSV2R1FKKb/xtWisF5GfAU4R6SIizwOL/ZiryekSk06ZQ1i48jO7oyillN/4WjTuAXoBZcAM4Chwn79CNUVDel4DwMqdc21OopRS/uPr1VMlwCPWQ9UgvctQor/3sLNMLy5TSjVfvnb32hV4AEiq+hpjzCj/xGp6HE4nya5Qsp3ndfuNSqlmzteb+z4EXsHbEZPbf3GatsSQZFaykV37tpCY0NXuOEopVe98PafhMsa8bIxZZoxZceLh12RNUGq7iwBYsOrDOuZUSqmmydei8W8RuVtE2opIzImHX5M1QaMyxuEwhk0Hl9odRSml/MLXw1MTrX9/W2WcQW/wO0lsVBsSK4Rdnhy7oyillF/4evVUsr+DNBcdiSMz6BAuVwUBAYF2x1FKqXp12sNTIvK7KsPjq037s79CNWWdI3tzzOFgydov7Y6ilFL1rq5zGjdWGX642rTLznal1jmRr0Vkq/VvdC3z/VVE1ovIRhF5TkTkbNfZUAZ2uwKAZVu0aCilmp+6iobUMlzT8zPxEDDPGNMFmGc9P3nhIkOBYUAa0BsYAAw/h3U2iIG9RhPh9rD96Aa7oyilVL2rq2iYWoZren4mrgXetobfBq6rZd0hQBAQDAQCB89hnQ0iICCQ5IoQdku+3VGUUqre1VU0+ojIUREpAtKs4RPPU89hva2NMSfaET8AtK4+gzFmCTAf2G895hhjNta0MBG5Q0QyRSTz8OHD5xCrfiQFdWR3oOFA7h67oyilVL06bdEwxjiNMS2NMRHGmABr+MTz014aJCJzRSSrhse11dZhqGGvRUQ6Az2A9kA7YJSIXFhLzleNMRnGmIz4+Pg63rL/9Ww7BCPCgpV6k59Sqnnx9ea+M2aMGW2M6V3D4zPgoIi0BbD+PVTDIsYAPxhjio0xxcCXwBB/5a1Pw/vdAEDWvu9tTqKUUvXLb0WjDp/z4w2DE4GaOqHYDQwXkQARCcR7ErzGw1ONTftWSXQoN+wq22V3FKWUqld2FY2pwCUishUYbT1HRDJE5HVrno+A7cA6YA2wxhjzbzvCno2OJpqdgcfxuLV9R6VU8+FrMyL1yhiTB1xcw/hM4DZr2A3c2cDR6k1KRA++L19C5ob5DEwdbXccpZSqF3btaTR7GZ0uBeCHTV/YnEQppeqPFg0/uSD9akI9HrYVrLU7ilJK1RtbDk+dD4KCgkmuCGI3uXZHUUqpeqN7Gn7UMaDdMIgYAAAe2ElEQVQdu4I8HCm0/4ZDpZSqD1o0/Kh7/EBcInyzQm/yU0o1D1o0/Ghk3/GIMWRmf2V3FKWUqhdaNPwouX0PepYHssKzQ+/XUEo1C1o0/GxA5FD2Bwr/+f7tumdWSqlGTouGn9006ncEeQxzN79ndxSllDpnWjT8LCE+kbSKMFY591NaVmJ3HKWUOidaNBrA0DY/4YjTwYfznrM7ilJKnRMtGg3gxot/Q4Tbw3c5Taa9RaWUqpEWjQYQERZFujuO1YEF5BUcsDuOUkqdNS0aDWRE8jiOOxzMmPeU3VGUUuqsadFoINePmEK8y8PS3AV2R1FKqbOmRaOBBAQE0peOZAWXsnPvJrvjKKXUWdGi0YCuSP0FLhFmLvir3VGUUuqsaNFoQKMGjKVjOWQeW2F3FKWUOitaNBqQiNA/qAebgz2s3PCt3XGUUuqMadFoYNcPvheAT5bpjX5KqaZHi0YDS+92Ad3LnKys2Kwt3yqlmhwtGjbICB/AniBh7rIP7I6ilFJnRIuGDW4c8VucxjBn/TS7oyilGpjbY/jz7I3cMT2TOesPUOH22B3pjATYHeB8lJjQlbTyFqxw5lBeXkZQULDdkZRSDaC0ws39M1fzZdYBokID+e+Gg8SFBzO2XzvGZ3Sgc6twuyPWSfc0bDI4biR5AQ4+WfCS3VGUUg2gqLSCSW8t48usA/zvVT3JfGQ0r9+SQd+OUby+aCejn/6WcS8v5oPMPRwrc9kdt1ZijLE7Q73KyMgwmZmZdseo05HCw1z2yQj6VcTw8h3f2R1HKeVHh4vKmPTWMjYfKOKp8X24rm+7k6YfKirlk5V7+WD5HnbkHiMsyMnVfRK4YUAH+naIQkT8nlFEVhhjMuqaTw9P2SQ6Mp70imhWBeZTWJxPZHiM3ZGUUn6wO6+En7+5lENHy3h9YgYjurU6ZZ5WESFMGd6JOy9KIXPXEWYu38Nnq/fx/vI9JMeFERsWVDlv9Z/5VX/4d20dwdSxaf56K4AenrLVRR2v4ZjDwftzteVbpRqrrL2FzN986KxOWK/fV8j1Ly+m8HgF790+qMaCUZWIMCAphqfG92H5H0Yz9fpUkuPCCA50VD5CAh20CHRWPkKDAggL9j5CAp1n+zZ9poenbFRaVsKl/xpAiiuMaXcuszuOUqqKknIXf/1qM28vycYYaBURzA0ZHfjpgA50iAmt8/VLtudxx/RMIkICmP6LgXRuFeHXvG5XOWWlBYSGn74w1UYPTzUBIcGh9DMJLAzez54DO+jQJsXuSEopYPG2XB78ZC178o8zcUgiQzvH8cHyPby0YBsvLtjGhV3i+dnADlzcozWBzlMP2HyVtZ97319Nx5hQpk8eSEJUizNav9tVTkFBNvkFO8g/uouCY4coLDnMkdIjFJQXUlheTIG7hEJ3GQXGRYEYigTSCWb6JP+2badFw2aXdv8587b/jZkL/soDN75idxylzmtFpRX8efYmZizbTVJsKB/cOYSByd7zjT/p1YZ9Bcf5IHMPHyzfw5R/rSQuPJjxGe25cUAHEmPDAJixbDePzFpHnw5RvDlxANHW+QiP20VBwU4O5W7kUMFOcov2kF9ymLzSI+SVF5LvLiHfU04ebgoEPLWc/A71GKKMEClOohxBtHO2JDIwjKigSDpGdfL7NtLDUzbzuN1c/lYfojyBzLxjld1xlDpvLdh8iIc/WcfBo6X84oJkfnNJN1oE1XyOwO0xLNxymPeW7eabTYdwmhIuTiqgTeh+Nu7dTOvIY7SOKSe3vJBDrhIOmXIOO8BVQyEI9RhijRDjCCTW2YKYwHBig6OJaRFLbFgbYsITiIpoR1TLDkRGdiQo2D+HufTwVBPhcDrp7+zCFwFbydq2lN6dB9kdSanzSmFJBU/8ZwMfrcihc6twPrprKP06Rp80T3HRPvYdWM2BvE3sK9jJvuK97C/NpcRdRGLnCnIdsEgEPEBbWA+ElRhaGQetnMH0D2hFq5BoWoW2plVEe+KjkomLSiE2tistWkTXmKux0qLRCFyX8Uv+vfJ+Pl78DL07v293HKXOG/9df4BHPs0i/1gpvxziZHhiNjnbvmLxip3sOXaAPRWF7KaCQsfJewiBxtDWI7R1hjAsKJaEFq1oHd6egMD29O6YRuv4HoSFt7bpXfmXFo1GYGDqaLosEZbJelyuCgICAu2OpFSztO/wLjI3zmXb/lVkF2zjsCuXNq3KcAR4mF7gYHqBdz6HVRQ6Olvwk+DWtA9vR9uWibSN7UZCfG9iY7vhcJ6fX5/n57tuhEZEj+K1knm8OOsBfj3+H3bHUarJ8rjdbN61mjXbvmXn4bXsL9nNIc8RDgSUkxfw45VOAcGGhAAhMTCcYSFxdIjoSMfYbnRs3ZeENv0IDA6z8V00Xlo0Gom7x/yN+W/2Z5ZrLjcf2U9cdFu7IynV6B3I3cMP62azaf8yco7tYL/JZ2+gi2OOH4tDaICHhAon3dwxtAlIoGNMT3omDaFv1wsICa77fgt1Mi0ajURAQCCTuv6KP+x6nr/PmsJfJn9mdySlGg2Xq4LMjfNZs30B2UfWs698H/ucxzkQ+OO5hnCnhw6uQAa4WtM2NJGU+DTSOl9E98R0HE7/3yl9vtCi0YhcO+IO/v3qW3wduJ2btiwmretQuyMp1eBcrgqWrZ/Lii3/ZUfhenI8h9gVWMFxa+/BIYYEJ3TwhDPYdCAlLo3+3UbTO2WgFocGYEvREJHxwGNAD2CgMabGGytE5DLgH4ATeN0YM7XBQtrk3hF/Z/L3d/Dc/P/h9a5L7I6jlF953G5WbVrIss1fsj1/HTmeA2QHllceXgp2GpI8Toa42pLYsis92g9hSNrlREXE2Zz8/GXXnkYWcD3wz9pmEBEn8CJwCZADLBeRz40xGxomoj3Sug7lkkUpfBGyk88WvMq1I+6wO5JS9eZA7h6+WfE+6/ctJrt8FzsDyyiymuEIDDAklTsY6GpNUsse9E0ZxZC0y/W8QyNjS9EwxmwE6mojfiCwzRizw5r3feBaoFkXDYDfjHmFJZ9cwttbXuDKC27VS3BVk+Rxu1m+4RuWbPw32wrXsos8dgUajAgiho4Ooa87jqQW3emTNIJhfa4iLNS/jfqpc9eYz2m0A/ZUeZ4D1Hi7tIjcAdwB0LFjR/8n87P46ASui7iYN47P10twVZNx7Pgx5i6bwcrsr9lZuoOdgSUUWHsRYQEeOpWHcJXpSO+2QxnV/0baxHWwObE6G34rGiIyF2hTw6RHjDH1emmQMeZV4FXwtj1Vn8u2y6+u/zvfWpfgTig4QGxUTZtSKfsUHSvg66Xvs3LX12yv2Mm2wHJKrTun2zkNvVwxdGrRk4FdL2dYnyt1j7mZ8FvRMMaMPsdF7AWq/hRpb407L1S9BPepT+7US3CV7QqL85mz5F+szpnPjopstgVWUOYQcECiwFBXa3rGD2J0/5vp1KGX3XGVnzTmw1PLgS4ikoy3WNwI/MzeSA1LL8Ft/jweg8cYAmrok8Fu5eVlfL10Bku2/5ut5TvYFlhBucN7PiJJhAvdCfSKHcwlGRNITOhqd1zVQOy65HYM8DwQD/xHRFYbY34iIgl4L629whjjEpFfAXPwXnL7pjFmvR157XTiEtznF/wPr+kluM2Kx2O4451M1u0t5IWf9WNAkr39xHvcbpatn8eCrJlsKlrHlsBj3iubBJIFhrvb0TtuGJcMmKAdhp3HtD+NJuDhN6/hC+dO/pR4L9eOuN3uOKqevLFoJ098sYGo0ECKS108fEUPJg9Lquuqwnq1Y08WXy57m6y8H9jizOeQ1TZTK5eHbu5YescN5rIBE0nRw03Nnq/9aWjRaAIOH9nH+E8uIcbt5IPJK/SEYjOwfl8hY15czEVd4/j7Dek88OEavt5wkCvT2vLk2DTCg/1zEKC0rIQ5S95hyY4v2OTaxfZg7///CLeHbhXhdIvozajUm8joOVLvrj7PaCdMzUjVS3BfmvUA9+oluE3a8XI3v35/NVGhgfx1XB8iWwTyzwn9+efCHfxtziY2HyjilQn96Nyqfu5Z2JK9mv8sfZ2sgkw2BhZR5HTgcBi64ORaTxeGdLmaSwbeRFBQcL2sTzVvuqfRRLhcFYx/sz9HnG4+HjtPL8Ftwh6ZtY53l+7mX78YxAVdTm4OY/H2XO6dsYrj5W6eHJfGVWkJZ7z80rIS/vP9NJbu/A+bPLvZ6e2imliXh+7uGPrEDeWqoXfqeQl1Et3TaGZOugR31p385Va9BLcpmrP+AO8u3c2dF6WcUjAAhnaK44t7LuTud1fwq/dWsXJXAQ9f0Z3AOq6u2rN/K59+/zJr8pewIeAoRU4HTqehqzuAMXTlwu7juDhjrB5yUudM9zSamNteHcLqwCL+OfgV+ve4wO446gwcKCzlsn8spH10Cz65axhBAbUXgnKXhz/P3si0xdlkJEbz4s39aN0ypHK6x+1m8dovmbf2X6wv28SWIBduEaLdHnq6Yujb6iKuHnYHCfGJDfHWVDOgJ8KbqbVbFnP797cT6xLu7T+Ny/rV+TdWjYDHY5jwxlJW7S7gP/deQEp8uE+v+3zNPh76eC2hQQH8/fou7NvzEUt3f8UG2c9+qy+JpHLo6Uzmgi5juHzIBL1QQp0VPTzVTKV1Hcr/Hbyb/932Mi8vv5W1+17mN5cPbZQ3h6kfvfrdDhZvz+PJsak+FwyA4UkuHhmygDnZ3/C7pccocTgICjD0KA9hVGg6Vwy4jbQug/2YXKmTadFogq648JcEOo7x263TCdjzK37+2l959mcXnXT4QjUea3MKeGrOZq5IbcMNGXU30rd371Lmr32LBQdXsILjuESIDTYMMfGkxYzkuovuIiayVQMkV+pUWjSaqEuG/Y4/VRTx++xZ9Cp5mKuffZynbxpc48lVZZ9jZS7unbGKVhHB/GVMWo037hm3mw1bZvHNpg9ZcGQTW5weAFKMMLFld0Z2H09q97E4nPrfVdlPP4VN2FUjnqD4qwL+38EF9OdP3PLmw9wzqgf3XtwFp6Ph7ipWtXvs8/Xsyi/h/dsHExn647mGirJilq2Zxvwd/2F+yR4OOQWHMfR1hPBAXF9G9L6FxMQLbUyuVM20aDRxN172PMWfT+QfrOTy5Gf5x7zfkLkrn2d/2pf4CL1Zy05frN3HhytyuGdUZwalxFJUmMN3q19j/p75LKrIp9ghtPAYhgZGMbLdhVzU93aio/XeCdW46dVTzYDxeHhm1jjeKt7KGElm5pa7aNkikOdu7MuQTrF2xzsv5Rwp4fJ/fEd63EEu67SUbw8tZ7kpwSVCjMcwMiSBkSmXMyjtVkJaRNkdVym95PZ8YzwenvjgCj4s28vtYel8uusXZOce44GfdOOu4Z0atBG885nH7WbBill8vPxtcgKz2WHt7CW5YWTLLozqNpbUHuNxBgTZG1SpavSS2/OMOBw8Mu5zit+/hNeOrebhXrNYcmQSf/1qMwcKS3ns6l449DyHX5SUHuPzhf/kh12zyZL9HAx0QDj0dDn5dUw6o3rfQkryKLtjKlUvtGg0I86AIP7fT7/k2IyLmXpgHlOTo2gbeSOvLtxBQUkFT43vc9q7kJXvcg5l89l3L7A693uyAo9S7HQQHGDoWRHKZcH9uWbIFLomptsdU6l6p0WjmQkMDOXv42dz18xLeGTnxzzTLYKYy8cw9ctNFB6v4OUJ/QgN0j/72Vix4Vu+WvkmWSVr2RRUgUuEqEAP6a5YMuJGMWb43Xr/hGr29JxGM1VctJ/bPrqcDeJiXEg7Uto+xqNf5ZPeIYo3Jw0gKlSPqdfl2PFivlj0Ost2fclGs5c9Qd7Dex3KDb0cHbmg8xguH3KLNimumgU9Ea44VnyAl768k3ePbaelgRsiRvL8qp+QGBfB9MmDaBOpd5BXt2XXar744TXWFWSyMaiYYw4HgcbQvSyInqG9uCT9FgalXmJ3TKXqnRYNVWnzli/40+JHWS0VpLkDyd83nqLAQbzzi0Ekx4XZHc9W5eVlzPnhXRZv/4yNFTsre7KLdXno5Ymjb+vhXHvhFOKjz7xfC6WaEi0a6iQet4vP5j/M03u+pEig39FWbD16Fy/feim920XaHa/BeNxulq7/mm/XfcCm4iy2BB7z9mRnDF3KnfQI6spFPcZr3xPqvKNFQ9Wo4MhOnv3qTj4u30+8y0PbvKH86vo/M6RzvN3R/GbLrtXMWT6d9fnL2eLM57B1BVkrl4du7lh6xQ3iqsG3k5jQ1eakStlHi4Y6rdVZ7/HEsifZ4vTQsySAG3s/wbXDriT3WBkHC8s4cLSUA4XHrX/LOHi0lANHSzlYWEpkaCB3XpTC+IwOhAQ2vl/jW3atZeGaj1h/6Ae2mv3sss75R7g9dHeF0y08jVF9bqR/9xG6N6GURYuGqpOropTpX/2af+YuwoXQvSQET1E3NhWN4ojb2we50yHEhwfTOjKENi2DadMyhKx9R1mx6witWwYzZXgnbhrY0bbiUVJ6jO9WzGJl9lx2HNtCtuMoB6zOiYI8hq4VQXQN7sKQLtcwKmOcXumkVC20aCif7cpZzTNz/oc1coBc6/h+b08gQ6N7Mzr1ZrqmXIo4frwp0BjDku15/GPeVpbuzCcuPJgpw1P42aCOfr0HxON2s23POhauncXmw8vZ5d7HjkAXZdad7rEuD8mucJJadKJPxxGMzBhHZHiM3/Io1Zxo0VBnzON2sXHr53y78UMWFGxko8MNQDs3DA/vyPCUKxmQeguBwT/2PPfDjjye/2Yr32/LIzYsiNsvSuHngxMJCz774lFaVsLaLd+Tlb2YPfkbOVC2l8McZX+Ai6NWD4UBxpBc7qCjozVdY/pyQe/r6N1pkB5uUuosadFQ5+zgwbUsXPMW3x74gR/cRZQ5hDCPIc0RSnxAODFBLYkJiSEmrBXHK6JZnO1k6Z4WeALbM+miXtwyJJGIEG8fEuXlZRwq2Efukb3kFx6goOgQR0tyKSo9wrHyQo6UHeJwxWEOOo6xPxAqqjSwGOPy0MYdRLxE0SakA73bD2NExjiiIrQFX6XqixYNVa+Ol+SzbO00FmTPYdPxQ+QbF3liKg8NVRfu9hDhESrEUOKAEsfp27xyGkNbF7Ryt6BVQDxtw1Po3LYv/bpfTPtWSX54R0qpqrRoKL8zHg/HS3LJO7KdvMJs8o/mkH9sP3klh9lbnMv+Y4U4cdLC0YIQZyihAeGEBbUkLDiKiBaxRIXFE92yFbGRCbRv3YnQkPP7RkOl7KRNoyu/E4eD0PBWhIa3okOHIXbHUUo1AG0nWymllM+0aCillPKZFg2llFI+06KhlFLKZ1o0lFJK+UyLhlJKKZ9p0VBKKeUzLRpKKaV81uzuCBeRw8Cuc1hEHJBbT3H8rSllhaaVtyllhaaVtyllhaaV91yyJhpj6uyNrdkVjXMlIpm+3ErfGDSlrNC08jalrNC08jalrNC08jZEVj08pZRSymdaNJRSSvlMi8apXrU7wBloSlmhaeVtSlmhaeVtSlmhaeX1e1Y9p6GUUspnuqehlFLKZ1o0lFJK+ey8KRoicpmIbBaRbSLyUA3TfyMiG0RkrYjME5HEKtPcIrLaenzeSPJOEpHDVXLdVmXaRBHZaj0mNoKsz1TJuUVECqpMa9BtKyJvisghEcmqZbqIyHPWe1krIv2qTGvQ7epj3putnOtEZLGI9KkyLdsav1pE/N6dpQ9ZR4hIYZW/96NVpp32M2RD1t9WyZllfU5jrGkNul2tdXYQkfnWd9R6Efl1DfM0zGfXGNPsH4AT2A6kAEHAGqBntXlGAqHW8F3AzCrTihth3knACzW8NgbYYf0bbQ1H25m12vz3AG/auG0vAvoBWbVMvwL4EhBgMLDUju16BnmHnsgBXH4ir/U8G4hrRNt2BPDFuX6GGiJrtXmvBr6xa7ta62wL9LOGI4AtNXwnNMhn93zZ0xgIbDPG7DDGlAPvA9dWncEYM98YU2I9/QFo38AZq6oz72n8BPjaGJNvjDkCfA1c5qeccOZZbwJm+DHPaRljFgL5p5nlWmC68foBiBKRtjT8dvUprzFmsZUHbP7c+rBta3Mun/ezcoZZbf3MAhhj9htjVlrDRcBGoF212Rrks3u+FI12wJ4qz3M4dYNX9Qu8FfuEEBHJFJEfROQ6fwSsxte8Y63d0I9EpMMZvra++Lw+65BfMvBNldENvW3rUtv7aejtejaqf24N8F8RWSEid9iUqbohIrJGRL4UkV7WuEa7bUUkFO8X7MdVRtu6XUUkCegLLK02qUE+uwFn+8LmSkQmABnA8CqjE40xe0UkBfhGRNYZY7bbk7DSv4EZxpgyEbkTeBsYZXOmutwIfGSMcVcZ1xi3bZMjIiPxFo0Lqoy+wNq2rYCvRWST9QvbLivx/r2LReQK4FOgi415fHE18L0xpupeiW3bVUTC8Raw+4wxRxtindWdL3sae4EOVZ63t8adRERGA48A1xhjyk6MN8bstf7dASzAW+X9qc68xpi8KhlfB/r7+tp6dibru5Fqu/k2bNu61PZ+Gnq7+kxE0vB+Bq41xuSdGF9l2x4CZuE9DGQbY8xRY0yxNTwbCBSROBrxtuX0n9kG3a4iEoi3YLxrjPmkhlka5rPbkCdz7Hrg3aPagffQyIkTbb2qzdMX78m4LtXGRwPB1nAcsBX/n6TzJW/bKsNjgB/Mjye9dlq5o63hGDuzWvN1x3sCUezctta6kqj9ZO2VnHwycZkd2/UM8nYEtgFDq40PAyKqDC8GLrM5a5sTf3+8X7S7re3s02eoIbNa0yPxnvcIawTbVYDpwLOnmadBPrvnxeEpY4xLRH4FzMF7pcabxpj1IvJHINMY8znwNyAc+FBEAHYbY64BegD/FBEP3j2zqcaYDY0g770icg3gwvvBnmS9Nl9EngCWW4v7ozl519qOrOD9xfa+sT7FlgbftiIyA+9VPHEikgP8HxBovZdXgNl4r0LZBpQAt1rTGnS7nkHeR4FY4CXrc+sy3lZOWwOzrHEBwHvGmK9szjoOuEtEXMBx4Ebr81DjZ8jmrOD9MfZfY8yxKi9t8O1qGQb8HFgnIqutcb/H+6OhQT+72oyIUkopn50v5zSUUkrVAy0aSimlfKZFQymllM+0aCillPKZFg2llFI+06KhGj0RKfZhnvusJh/qa53XiUjPelze4nN4bbH1b4KIfHSa+aJE5O6zXY9SvtCioZqL+4AzKhoi4jzN5OuAeisaxpih9bCMfcaYcaeZJQrQoqH8SouGajKs/hgWWA00bhKRd60+BO4FEoD5IjLfmvdSEVkiIitF5EOrzZ4TfSE8KSIrgfEicruILLca0ftYREJFZChwDfA3q8+ETiKSbjWquFZEZolItLW8BeLtLyRTRDaKyAAR+cTqt+BPVbIXVxl+ULz9MawRkak1vM9kK/u6astIEqv/BxHpJSLLrHxrRaQLMBXoZI37m4iEi7dvmJXWsq6tspyNIvKaePtm+K+ItLCmdRaRuVa2lSLSyRr/W2s7rRWRx+v1D6uaFn/f/q4PfZzrA6vPDbx38BbibTvHASzB23gcVOnjAG+TJAuxmn8AHgQerTLf76osO7bK8J+Ae6zhacC4KtPWAsOt4T9iNeeAt72sJ63hXwP78PZ9EIy3NdHYau/hcrxNT5zou+WU5hyAz4FbrOFfVnltElazF8DzwM3WcBDQgmrNYuC9Y7lllW2yDW8TE0l4WxJIt6Z9AEywhpcCY6zhELx7b5cCr1qvdQBfABfZ/bnQhz2P86IZEdWsLDPG5ABYzSkkAYuqzTMY76Gl763mHoLwFpgTZlYZ7m39mo/C24zMnOorFJFIIMoY86016m3gwyqznGgqZR2w3hiz33rdDrwNxeVVmXc08Jax+m4xNTfnMAwYaw2/AzxZwzxLgEdEpD3wiTFmq/VeT4oO/FlELgI8eJvDbm1N22mMOdEcxQogSUQigHbGmFlWtlLrfVyKt3CssuYPx9s6rZ0t5iqbaNFQTU1ZlWE3NX+GBW+nMzfVsoyqbQlNA64zxqwRkUl492bONpOnWj5PLfl8cdr2fYwx74nIUryN1M0Wb/P4O6rNdjMQD/Q3xlSISDbevYeqmcG7HVucZnUC/MUY888zyK+aKT2noZqLIrzdYIK3B7thItIZQETCRKRrLa+LAPaLt9npm2tanjGmEDgiIhda034OfMvZ+Rq49cSVXmL1O13N93gbeKRapkri7X9khzHmOeAzII2TtwF4W2k9ZBWMkUDiqUv6kfH2CJcjVmdYIhJs5ZwDTK5yXqidePuSUOchLRqquXgV+EpE5htjDuNt9XeGiKzFeyiney2v+1+8x/G/BzZVGf8+8FsRWWWdDJ6I98T4WiAd73mNM2a8LaJ+DmRah9ceqGG2XwO/FJF11N7D2g1AlrWM3ni7+czDe0guS0T+BrwLZFjLuaXa+6vNz/G2oLwW77mXNsaY/wLvAUusZX3EycVJnUe0lVullFI+0z0NpZRSPtOioZRSymdaNJRSSvlMi4ZSSimfadFQSinlMy0aSimlfKZFQymllM/+P5lf+sE2iTasAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fee141210f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "for value, bootstrap in conditions.items():\n",
    "    plt.plot(results[f\"points_{bootstrap}\"], results[f\"energies_{bootstrap}\"], label=f\"{bootstrap}\")\n",
    "plt.plot(results[\"points_np\"], results[\"energies_np\"], label=\"numpy\")\n",
    "plt.legend()\n",
    "plt.title(\"Dissociation profile\")\n",
    "plt.xlabel(\"Interatomic distance\")\n",
    "plt.ylabel(\"Energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare number of evaluations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "21482\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "3908\n",
      "np\n",
      "Total evaluations for np:\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "for condition, result_full in results_full.items():\n",
    "    print(condition)\n",
    "    print(\"Total evaluations for \" + condition + \":\")\n",
    "    sum = 0\n",
    "    for key in result_full:\n",
    "        if condition != \"np\":\n",
    "            sum += result_full[key].raw_result.cost_function_evals\n",
    "        else:\n",
    "            sum = 0\n",
    "    print(sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extrapolation \n",
    "\n",
    "Here, an extrapolator added that will try to fit each (param,point) set to some specified function and suggest an initial parameter set for the next point (degree of freedom). \n",
    "\n",
    "- Extrapolator is based on an external extrapolator that sets the 'window', that is, the number of previous points to use for extrapolation, while the internal extrapolator proceeds with the actual extrapolation.\n",
    "- In practice, the user sets the window by specifying an integer value to num_bootstrap - which is also the number of previous points to use for bootstrapping. Additionally, the external extrapolator defines the space within how to extrapolate - here PCA, clustering and the standard window approach. \n",
    "\n",
    "In practice, if no extrapolator is defined and bootstrapping is True, then all previous points will be bootstrapped. If an extrapolator list is defined and no points are specified for bootstrapping, then the extrapolation will be done based on all previous points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Window Extrapolator: An extrapolator which wraps another extrapolator, limiting the internal extrapolator's ground truth parameter set to a fixed window size\n",
    "2. PCA Extrapolator: A wrapper extrapolator which reduces the points' dimensionality with PCA, performs extrapolation in the transformed pca space, and untransforms the results before returning.\n",
    "3. Sieve Extrapolator: A wrapper extrapolator which performs an extrapolation, then clusters the extrapolated parameter values into two large and small clusters, and sets the small clusters' parameters to zero.\n",
    "4. Polynomial Extrapolator: An extrapolator based on fitting each parameter to a polynomial function of a user-specified degree.\n",
    "5. Differential Extrapolator: An extrapolator based on treating each param set as a point in space, and performing regression to predict the param set for the next point. A user-specified degree also adds derivatives to the values in the point vectors which serve as features in the training data for the linear regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Here we test two different extrapolation techniques"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define different extrapolators\n",
    "degree = 1\n",
    "extrap_poly = Extrapolator.factory(\"poly\", degree=degree)\n",
    "extrap_diff = Extrapolator.factory(\"diff_model\", degree=degree)\n",
    "extrapolators = {\"poly\": extrap_poly, \"diff_model\": extrap_diff}\n",
    "\n",
    "for key in extrapolators:\n",
    "    extrap_internal = extrapolators[key]\n",
    "    extrap = Extrapolator.factory(\"window\", extrapolator=extrap_internal)\n",
    "    # define extrapolator\n",
    "    # BOPES sampler\n",
    "    bs_extr = BOPESSampler(gss=me_gsc, bootstrap=True, num_bootstrap=None, extrapolator=extrap)\n",
    "    # execute\n",
    "    res = bs_extr.sample(es_problem, points)\n",
    "\n",
    "    results_full[f\"{key}\"] = res.raw_results\n",
    "    results[f\"points_{key}\"] = res.points\n",
    "    results[f\"energies_{key}\"] = res.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poly\n",
      "diff_model\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Energy')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4VGX6//H3PSWNhCRA6IGAdAgECEUjTbALqAi6wkrEhrvrKvvTr7ruurquLq6u69orIoqo6FoXdQVBQVAM0nsxVClJSG9Tnt8fM4khBDJAJiflfl3XXMyc+pmTYe45zznnOWKMQSmllAqEzeoASiml6g8tGkoppQKmRUMppVTAtGgopZQKmBYNpZRSAdOioZRSKmBaNFS9IiIviMifLVjvZyIy9TTn7SAi+SJir+lcNU1ErhCRvf68/UVko4iM9I97QETetDiispjodRqqrhCRdKAV4AY8wCZgDvCSMcZrYbRT4n8fNxpjFlqd5VSJyE7gD8aYj6oY9wDQxRgzpdaDqTpD9zRUXTPWGBMFdARmAncDr1obqWEQEUcAk3UENgY7i6q/tGioOskYk2OM+Ri4GpgqIn0ARGS2iPzN/7yFiHwqItkikiUiS0XE5h93t4jsF5E8EdkqIqP9w0NF5EkROeB/PCkioWXrFZHxIrJGRHJFZKeIXOQfvkREbvQ/P0tEvhKRTBHJEJG5IhLjH/cG0AH4xN/E838ikiAipuxLW0TaisjH/sw7ROSmCut/QETeFZE5/uwbRST5RNvJv9zfi8guf5bHKmyDVBH5VkT+JSKZwAMiYhORP4nIbhE57F9PtH+75AN2YK1/jwMRSReRMSdY91ARWe7f/mvLmrFUw6ZFQ9VpxpiVwD5gWBWj/59/XBy+Zq0/AkZEugO/Awb591ouBNL989wHDAWSgH7AYOBPACIyGF9z2F1ADDC8wnwVCfB3oC3QE4gHHvDn/TWwB98eU6Qx5h9VzP+2P3db4CrgERE5r8L4cf5pYoCPgWeq3jrlrgCSgQHAeGBahXFDgF34ts/DQKr/MQroDEQCzxhjSowxkf55+hljzjrZCkWkHfBf4G9AM+BO4H0Riasmq6rntGio+uAAvi+mylxAG6CjMcZljFlqfAfpPEAo0EtEnMaYdGPMTv88k4G/GmMOG2OOAA8Cv/aPuwGYZYz50hjjNcbsN8ZsqbxSY8wO/zQl/mU8AYwI5I2ISDyQAtxtjCk2xqwBXgGuqzDZMmPMAmOMB3gDX3E7mUeNMVnGmD3Ak8CvKow7YIx52hjjNsYU+d//E8aYXcaYfOBe4JoAm64qmgIs8Of0GmO+BNKAS05xOaqe0aKh6oN2QFYVwx8DdgD/8zfP3AO+L3XgDny//g+LyNsi0tY/T1tgd4Vl7PYPA98ew06qISKt/MvcLyK5wJtAiwDfS1sgyxiTVylDuwqvD1Z4XgiEVfOlvrfSstqeYFzZ+iu/fwe+PZFT0RGY6G+ayhaRbOBcfEVcNWBaNFSdJiKD8H2hLqs8zhiTZ4z5f8aYzviadP5QduzCGPOWMeZcfF9uBnjUP9sB/7AyHfzDwPcFe9JmGb9H/MtMNMY0xferWypGO8m8B4BmIhJVKcP+ANZ7IvGVlnWgwuvKWap6/27g0Cmucy/whjEmpsKjiTFm5ikuR9UzWjRUnSQiTUXkMnxt+28aY9ZXMc1lItJFRATIwdcs5RWR7iJynv8AdzFQBJSdsjsP+JOIxIlIC+B+fHsK4DtL63oRGe0/YNxORHpUES8KyAdy/G37d1Uafwjf8YLjGGP2AsuBv4tImIj0xdcsdibXP9wlIrH+pq/bgXdOMu08YIaIdBKRSHwF8B1jjPsU1/kmMFZELhQRu/+9jBSR9qf3FlR9oUVD1TWfiEgevl+y9+E7XnD9CabtCizE9wW+AnjOGLMY3/GMmUAGvqaelvja7sF34DYNWAesB370Dys76H498C98Rehrjv1VXuZBfAedc/AdDP5PpfF/x1eYskXkzirm/xWQgO9X/wfAX87wmo6PgFXAGn+ek52iPAvfcZJvgJ/wFdXbTnWF/uI3Ht/JB0fw/b3uQr9TGjy9uE+pekxEDNDVfxxHqaDTXwVKKaUCpkVDKaVUwLR5SimlVMB0T0MppVTATvUq0DqvRYsWJiEhweoYSilVr6xatSrDGFNtNzANrmgkJCSQlpZmdQyllKpXRGR39VNp85RSSqlToEVDKaVUwLRoKKWUCliDO6ahlDozLpeLffv2UVxcbHUUFQRhYWG0b98ep9N5WvNr0VBKHWPfvn1ERUWRkJCAry9I1VAYY8jMzGTfvn106tTptJahzVNKqWMUFxfTvHlzLRgNkIjQvHnzM9qL1KKhlDqOFoyG60z/tlo0/HZtXMesm/7EJ889bXUUpZSqs7Ro+NnsNors55Gx7mD1Eyul6qzU1FTee++9M17OmjVrWLBgwUmnSU9P56233jrjdQXihRdeYM6cObWyrpPRouGX0KMPdlc+piTc6ihKqTrgTIuG232qN0M8uenTp3PdddfV6DJPhxaNChzuLIwnxuoYSjVq6enp9OzZk5tuuonevXtzwQUXUFRUBPi+yIcOHUrfvn254oorOHr0aJXLWLhwIcnJyXTr1o1PP/0U8B3gv/7660lMTKR///4sXrz4hMNLS0u5//77eeedd0hKSuKdd97h66+/JikpiaSkJPr3709eXh733HMPS5cuJSkpiX/961/Mnj2bcePGcd555zF69Gjy8/MZPXo0AwYMIDExkY8++qj8Pfbo0YPJkyfTs2dPrrrqKgoLCwFfV0j/93//R2JiIoMHD2bHDt/9tR544AEef/xxAEaOHMndd9/N4MGD6datG0uXLgWgsLCQSZMm0atXL6644gqGDBlS490q6Sm3FYjJxGtva3UMpeqMBz/ZyKYDuTW6zF5tm/KXsb1POs327duZN28eL7/8MpMmTeL9999nypQpXHfddTz99NOMGDGC+++/nwcffJAnn3zyuPnT09NZuXIlO3fuZNSoUezYsYNnn30WEWH9+vVs2bKFCy64gG3btp1w+F//+lfS0tJ45plnABg7dizPPvssKSkp5OfnExYWxsyZM3n88cfLC9Ps2bP58ccfWbduHc2aNcPtdvPBBx/QtGlTMjIyGDp0KOPGjQNg69atvPrqq6SkpDBt2jSee+457rzTd3fg6Oho1q9fz5w5c7jjjjvKl1+R2+1m5cqVLFiwgAcffJCFCxfy3HPPERsby6ZNm9iwYQNJSUln9Leqiu5pVCDOXFzOZhQVFFgdRalGrVOnTuVfeAMHDiQ9PZ2cnByys7MZMWIEAFOnTuWbb76pcv5JkyZhs9no2rUrnTt3ZsuWLSxbtowpU6YA0KNHDzp27Mi2bdtOOLyylJQU/vCHP/DUU0+RnZ2Nw1H1b+7zzz+fZs2aAb7rIv74xz/St29fxowZw/79+zl06BAA8fHxpKSkADBlyhSWLVtWvoxf/epX5f+uWLGiyvVceeWVx2wfgGXLlnHNNdcA0KdPH/r27VvlvGdC9zQqsDUpxZQ62bD8awadf4nVcZSyXHV7BMESGhpa/txut5c3TwWq8mmlNXEK8T333MOll17KggULSElJ4YsvvqhyuiZNmpQ/nzt3LkeOHGHVqlU4nU4SEhLKr5E4WcYTPa+obBvZ7fYaP35yMrqnUUFEK98fe8/aDRYnUUpVFh0dTWxsbHn7/RtvvFG+11HZ/Pnz8Xq97Ny5k127dtG9e3eGDRvG3LlzAdi2bRt79uw56fCoqCjy8vLKl7lz504SExO5++67GTRoEFu2bDlumspycnJo2bIlTqeTxYsXs3v3L72P79mzp3wv4q233uLcc88tH/fOO++U/3v22WcHvI1SUlJ49913Adi0aRPr168PeN5A6Z5GBW16duPIXsjfn211FKVUFV5//XWmT59OYWEhnTt35rXXXqtyug4dOjB48GByc3N54YUXCAsL4ze/+Q233noriYmJOBwOZs+eTWho6AmHjxo1ipkzZ5KUlMS9997LsmXLWLx4MTabjd69e3PxxRdjs9mw2+3069eP1NRUYmNjj8kxefJkxo4dS2JiIsnJyfTo0aN8XPfu3Xn22WeZNm0avXr14tZbby0fd/ToUfr27UtoaCjz5s0LePv85je/YerUqfTq1YsePXrQu3dvoqOjT3Ern1yDu0d4cnKyOd2zBXIyM3jzj2sI937NtJcfquFkStUPmzdvpmfPnlbHaNDS09O57LLL2LDh+FaNshvJtWjR4pSX6/F4cLlchIWFsXPnTsaMGcPWrVsJCQk5Zrqq/sYissoYk1zdOnRPo4Lo5i1wlmaDN8rqKEopdcoKCwsZNWoULpcLYwzPPffccQXjTGnRqMTuzcQQW/2ESil1mhISEqrcywDKz4Q6HVFRUUG/3bUeCD9OFh57c6tDKKVUnaRFoxIJKcAVEsORA3utjqKUUnWOFo1KnNG+EwPWf73Y4iRKKVX3aNGoJKq9r2nq8LZ0a4MopVQdpEWjks6D+gNQfETvj6yUVdLT0+nTp88ZL+fJJ58s7wjwRGbPns2BAwfOeF2BOOecc2plPcGkRaOSXoPPxeYpxVMYWv3ESqk67UyLhsfjqdE8y5cvr9HlWUGLRiUOpxOnKxM8Ta2OolSj5na7j+s6fNGiRfTv35/ExESmTZtGSUkJQJXDn3rqKQ4cOMCoUaMYNWoUHo+H1NRU+vTpQ2JiIv/617947733SEtLY/LkySQlJVFUVERCQgJ33303AwYMYP78+bz88ssMGjSIfv36MWHChPIilJqayvTp04/rgn327NmMHz+ekSNH0rVrVx588MHy9xQZGQnAkiVLGDlyJFdddVV5F+llF1ovWLCAHj16MHDgQH7/+99z2WWX1eZmr5Yl12mISDPgHSABSAcmGWOq7BhfRJoCm4APjTG/q5V83kyM6Gm3SvHZPXCwhvsvap0IF8+sdrLKXYc/8cQTvPjiiyxatIhu3bpx3XXX8fzzzzN9+nRSU1OPG37HHXfwxBNPsHjxYlq0aMGqVavYv39/+fUR2dnZxMTE8Mwzz/D444+TnPzLxdDNmzfnxx9/BCAzM5ObbroJgD/96U+8+uqr3HbbbUDVXbADrFy5kg0bNhAREcGgQYO49NJLj1k+wOrVq9m4cSNt27YlJSWFb7/9luTkZG655Ra++eYbOnXqVN7bbV1i1Z7GPcAiY0xXYJH/9Yk8BFTd/3Gw2HNwO5vjdrlqdbVKqV9U7jp80aJFdOrUiW7dugG/dI2+devWKodX1rlzZ3bt2sVtt93G559/TtOmJ25NuPrqq8ufb9iwgWHDhpGYmMjcuXPZuHFj+biqumAHX/fozZs3Jzw8nCuvvPKYbs/LDB48mPbt22Oz2UhKSiI9PZ0tW7bQuXNnOnXqBFAni4ZVV4SPB0b6n78OLAHurjyRiAwEWgGfA9X2iVJTbOFFeLzh7Fq/mm4DBtfWapWqewLYIwiWyl2Cx8TEkJmZedrLi42NZe3atXzxxRe88MILvPvuu8yaNavKaSt2b56amsqHH35Iv379mD17NkuWLDlhxrLXgXTNXrn799rs3vxMWLWn0coY87P/+UF8heEYImID/gncWd3CRORmEUkTkbQjR46ccbjQZk4Atq1cecbLUkqdnspdhycnJ5Oenl7eBFTWNXr37t2rHA4c03V5RkYGXq+XCRMm8Le//a28+am67s3z8vJo06YNLpervAv1MlV1wQ7w5ZdfkpWVRVFRER9++GH5HlN1unfvzq5du8q7EinrIr0uCdqehogsBFpXMeq+ii+MMUZEqupq9zfAAmPMvupuoGKMeQl4CXy93J5e4l80P6s9RzMg+6dDZ7oopdRpqtx1+FNPPcXQoUOZOHEibrebQYMGMX36dEJDQ3nttdeOGw5w8803c9FFF9G2bVuefPJJrr/+erxeLwB///vfgV8OaIeHh1d5l7yHHnqIIUOGEBcXx5AhQ44pMFV1wQ6+pqcJEyawb98+pkyZctzxjBMJDw/nueee46KLLqJJkyYMGjTojLZhMFjSNbqIbAVGGmN+FpE2wBJjTPdK08wFhgFeIBIIAZ4zxpzs+McZdY1eZu/2rXz8z/2Es4hpLzx8RstSqr7RrtEDk5qaymWXXcZVV111zPDZs2cfc2/xU5Wfn09kZCTGGH7729/StWtXZsyYURORy51J1+hWNU99DEz1P58KfFR5AmPMZGNMB2NMAr4mqjnVFYyaEt+1Ow5XHqY4ojZWp5RS5V5++WWSkpLo3bs3OTk53HLLLVZHOoZVB8JnAu+KyA3AbmASgIgkA9ONMTdalKuc3Z0JJsbqGEqpOmr27NlVDk9NTSU1NfW0lztjxowa37OoSZYUDWNMJjC6iuFpwHEFwxgzG5gd9GAVCFl47fG1uUqllKrz9IrwE3Hk4XI2o6igwOokSilVZ2jROAFHpAtjs7P2m0VWR1FKqTpDi8YJRLTy9RGzb/0mi5MopVTdoUXjBNr29p0BXPBzrsVJlFKq7tCicQJ9R4wG48WTZ9UJZkopVfdo0TiBqJgYQkqzMK5Iq6Mo1eikp6fTs2dPbrrpJnr37s0FF1xAUVERI0eOpOzi3YyMDBISEgDf6a+XX345559/PgkJCTzzzDM88cQT9O/fn6FDh5KVlQXAyJEjuf3220lKSqJPnz6sXLkSr9dL165dKeuCyOv10qVLF2qiS6KGSH9Gn4TNmwU0szqGUpZ5dOWjbMnaUqPL7NGsB3cPPq5/0uNs376defPm8fLLLzNp0iTef//9k06/YcMGVq9eTXFxMV26dOHRRx9l9erVzJgxgzlz5nDHHXcAUFhYyJo1a/jmm2+YNm0aGzZsYMqUKcydO5c77riDhQsX0q9fP+Li4mrk/TY0uqdxMnIUt13vq6GUFTp16kRSUhIAAwcOLO/E70RGjRpFVFQUcXFxREdHM3bsWAASExOPmbesu/Hhw4eTm5tLdnY206ZNY86cOQDMmjWL66+/vubfUAOhexonISH5uG3RHN67m5bxHa2Oo1StC2SPIFgqdx1eVFSEw+Eo73CwuLj4hNPbbLby1zab7Zhux6vqtjw+Pp5WrVrx1VdfsXLlyuN6s1W/0D2Nk3DG+D5c679ZbHESpRRAQkICq1atAuC99947rWWUdTe+bNkyoqOjiY6OBuDGG29kypQpTJw4EbvdXjOBGyAtGifRNN7XNHV4226LkyilAO68806ef/55+vfvT0ZGxmktIywsjP79+zN9+nReffXV8uHjxo0jPz9fm6aqYUnX6MFUE12jl9mwYhlfv15KhGMh1z/zSI0sU6m6riF3jT5y5Mjj7gdeJi0tjRkzZrB06VILktWu+tg1er3QI3kINk8JpjC0+omVUvXWzJkzmTBhQvmNmdSJ6YHwk3A4nThdmRhvtNVRlFI1oOL9vSu65557uOeeWrldT72nexrVEG8mxqan3SqlFGjRqJbYc3A5m+N2uayOopRSltOiUQ0JL8ZrD2PHuh+tjqKUUpbTolGNsBYhAGz/7juLkyillPW0aFSjRRffLV9z9mjnZUrVVRU7MlTBpUWjGr2HDQfAld2wrmdRSqnToUWjGm07dcXhysWURFgdRalGIz09nR49ejB58mR69uzJVVddRWFhIYsWLaJ///4kJiYybdo0SkpKjplv1qxZ5b3ZArz88svMmDGjtuM3aHqdRgDs7kyM0S7SVeNz8JFHKNlcs12jh/bsQes//rHa6bZu3cqrr75KSkoK06ZN44knnuDFF19k0aJFdOvWjeuuu47nn3/+mCIxadIkHn74YR577DGcTievvfYaL774Yo3mb+x0TyMAQhZeuxYNpWpTfHw8KSkpAEyZMoVFixbRqVMnunXrBsDUqVP55ptvjpknMjKS8847j08//ZQtW7bgcrlITEys9ewNme5pBECceZTam1GQm0OTpnp1uGo8AtkjCJbKXZjHxMSQmZlZ7Xw33ngjjzzyCD169NDOB4NA9zQCYI90g9hZt1S7SFeqtuzZs4cVK1YA8NZbb5GcnEx6ejo7duwA4I033mDEiBHHzTdkyBD27t3LW2+9VX7DJVVztGgEIKJ1FAD71m+2OIlSjUf37t159tln6dmzJ0ePHmXGjBm89tprTJw4kcTERGw2G9OnT69y3kmTJpGSkkJsbGwtp274tHkqAO0Te3I4HQoP5lkdRalGw+Fw8Oabbx4zbPTo0axevfq4aSt3RLhs2TI9aypIdE8jAH2HjQLjwZOvNVapuiw7O5tu3boRHh7O6NGjrY7TIOm3YACaNI0mpDQL442yOopSjUJCQgIbNmw45fliYmLYtm1bEBKpMrqnESCbJwuDnnarlGrctGgESGxZeBx6Xw2lVOOmRSNAElqI29mUg7t3Wh1FKaUso0UjQM4Y34VG679eYm0QpZSykBaNAEV3iAMgY+c+i5MopZR1tGgEqMugZACKj5RUM6VSqqY98MADPP7449x///0sXLgQgKVLl9K7d2+SkpIoKirirrvuonfv3tx1111BzzN79mx+97vfnfE09ZGechugrv0Hs8TzJaYkzOooSjVaf/3rX8ufz507l3vvvZcpU6YA8NJLL5GVlYXdbrcqXqOgRSNADqcTpysT49UOC1XjsfTdbWTsza/RZbaIj2TYpG7VTvfwww/z+uuv07JlS+Lj4xk4cCCpqalcdtllZGdn8+677/LFF1/w2WefkZeXR35+PgMHDuTee+/l6quvPm55qamphIeHs3r1ag4fPsysWbOYM2cOK1asYMiQIcyePRuAefPm8cgjj2CM4dJLL+XRRx8F4LXXXuPvf/87MTEx9OvXj9DQUACOHDnC9OnT2bNnDwBPPvlkee+8DZEWjVMg3kyMTU+7VSrYVq1axdtvv82aNWtwu90MGDCAgQMHlo+/8cYbWbZsGZdddhlXXXUV4OsWfc2aNSdd7tGjR1mxYgUff/wx48aN49tvv+WVV15h0KBBrFmzhpYtW3L33XezatUqYmNjueCCC/jwww8ZMmQIf/nLX1i1ahXR0dGMGjWK/v37A3D77bczY8YMzj33XPbs2cOFF17I5s0Nt586LRqnQBw5lNq743a5cDidVsdRKugC2SMIhqVLl3LFFVcQEeG7Y+a4ceNqZLljx45FREhMTKRVq1bl99ro3bs36enp7N69m5EjRxIX5zvxZfLkyeX37Kg4/Oqrry6/8nzhwoVs2rSpfB25ubnk59fs3lldYknREJFmwDtAApAOTDLGHK1iug7AK0A8YIBLjDHptRa0cp6IErzuULakfU+fs8+1KoZS6jSVNSnZbLby52Wv3W43ztP4Mej1evnuu+8IC2scxzutOnvqHmCRMaYrsMj/uipzgMeMMT2BwcDhWspXpbDmvg/ZzrQfrIyhVIM3fPhwPvzwQ4qKisjLy+OTTz6plfUOHjyYr7/+moyMDDweD/PmzWPEiBEMGTKEr7/+mszMTFwuF/Pnzy+f54ILLuDpp58uf11dE1l9Z1Xz1HhgpP/568AS4O6KE4hIL8BhjPkSwBhj+f5ey24dyToEuXurv3uYUur0DRgwgKuvvpp+/frRsmVLBg0aVCvrbdOmDTNnzmTUqFHlB8LHjx8P+E77Pfvss4mJiSEpKal8nqeeeorf/va39O3bF7fbzfDhw3nhhRdqJa8VxBhT+ysVyTbGxPifC3C07HWFaS4HbgRKgU7AQuAeY4yniuXdDNwM0KFDh4G7d+8OSu7De3cz/+GdhHsXMe2lh4OyDqWstnnzZnr27Gl1DBVEVf2NRWSVMSa5unmD1jwlIgtFZEMVj/EVpzO+qlVV5XIAw4A7gUFAZyC1qnUZY14yxiQbY5LLDlQFQ8v4jjhKczClkUFbh1JK1WVBa54yxow50TgROSQibYwxP4tIG6o+VrEPWGOM2eWf50NgKPBqUAIHyOHJBPQWkkrVVQ8//PAxxxwAJk6cyH333WdRoobFqmMaHwNTgZn+fz+qYpofgBgRiTPGHAHOA9JqL+KJZOG1JVgdQqmgMsbgazmuf+677z4tECdxpockrDp7aiZwvohsB8b4XyMiySLyCoD/2MWdwCIRWQ8I8LJFecuJM5/SkFjysrOtjqJUUISFhZGZmXnGXy6q7jHGkJmZeUanB1uyp2GMyQSOu4GvMSYN38HvstdfAn1rMVq17FFuKLKz7utFpIyfYHUcpWpc+/bt2bdvH0eOHLE6igqCsLAw2rdvf9rz6xXhp6hJm6bk74IDm7b6ThxWqoFxOp106tTJ6hiqjtKu0U9Rh6Q+ABQetPyyEaWUqnVaNE5RYsooxOvGk687aUqpxke/+U5ReJMmOF1ZGG9Tq6MopVSt0z2N02DzZGJoZnUMpZSqdVo0ToPYj+J2tsTtclkdRSmlapUWjdNgjy7E44hg+SfvWR1FKaVqlRaN09Au2XdjmvTlGyxOopRStUuLxmkYftU12N1FuDPCrY6ilFK1KqCiISL/EZFLRUSLDBASFoazdDdec/pXVSqlVH0UaBF4DrgW2C4iM0WkexAz1QsScpCSsHYc3L3T6ihKKVVrAioaxpiFxpjJwAB89/ReKCLLReR6ETn1m+o2ABHt7SB2vp3/vtVRlFKq1gTc3CQizfHdBOlGYDXwb3xF5MugJKvjksZeBEDujhyLkyilVO0J9JjGB8BSIAIYa4wZZ4x5xxhzG9Aob2PXY+AQQkoO4S1qYXUUpZSqNYF2I/KUMWZxVSMCuadsQ2Xz7sFt74bb5cLhbJStdEqpRibQohErIldWGpYDrDfGVHWr1kbB3jSHYlc0aQs/Y+jF46yOo5RSQRfoMY0bgFeAyf7Hy8DdwLci8usgZavz4hJ9p9zuWLLS4iRKKVU7Ai0aTqCnMWaCMWYC0AswwBB8xaNRGj7pWmyeUkqPaNOUUqpxCLRotDfGHKrw+jAQb4zJAhptr31RMTGElO7B69GL/JRSjUOgxzSWiMinwHz/6wn+YU2A7KAkqyfEvp8SxzkcPXKQ2LjWVsdRSqmgCnRP47fAa0CS/zEH+K0xpsAYMypY4eqD0DZgbE6+eXue1VGUUiroqt3TEBE7sNBfHPTy50q/Lwd/AAAgAElEQVR6nz+Mb9+GrM1HrI6ilFJBV+2ehjHGA3hFJLoW8tQ7SSPH4CzJwlugd/JTSjV8gR7TyAfWi8iXQEHZQGPM74OSqp6xe3fjsXW0OoZSSgVdoEXjP/6HqoKtSSbF3v6sXfoV/YadZ3UcpZQKmoCKhjHmdREJBzoYY7YGOVO9E9s9jsLNsPGLr7VoKKUatEA7LBwLrAE+979OEpGPgxmsPhl29STE66b4Z6uTKKVUcAV6yu0DwGD812QYY9YAnYOUqd5p3rodISV7Ma62VkdRSqmgCrRouIwxlW8c4a3pMPWZzbaf0tCOFOTq/TWUUg1XoEVjo4hcC9hFpKuIPA0sD2KuesfZshSvPYSv33nL6ihKKRU0gRaN24DeQAkwD8gF7ghWqPqoy3DfbUUOr9trcRKllAqeQO8RXmiMuc8YM8gYk+x/XhzscPXJoAsuxVGaiydXr4FUSjVcAZ1yKyLdgDuBhIrzGGP0/FI/h9OJw52O197B6ihKKRU0gV7cNx94Ad+NmDzBi1O/ScQRSunLth9X0m3AYKvjKKVUjQu0aLiNMc8HNUkDEH1WNEU74cdPPteioZRqkAI9EP6JiPxGRNqISLOyR1CT1UMpkyaA8VK4t9Hel0op1cAFuqcx1f/vXRWGGfQCv2O07ngWocXLMEZvxqSUapgC7XuqU7CDNBQ2215K7f0pLS4mJCzM6jhKKVWjTto8JSL/V+H5xErjHglWqPrMEVuIxxHO0vfftjqKUkrVuOqOaVxT4fm9lcZddLor9R8T+VJEtvv/jT3BdP8QkY0isllEnhIROd111pYOQ3sDsO+HbRYnUUqpmldd0ZATPK/q9am4B1hkjOkKLPK/PnbhIucAKUBfoA8wCBhxBuusFeeMm4DdXYAnO9LqKEopVeOqKxrmBM+ren0qxgOv+5+/Dlx+gnWHASFAKOAEDp3BOmtFSFgYztLdeIm3OopSStW46opGPxHJFZE8oK//ednrxDNYbytjTNndJw4CrSpPYIxZASwGfvY/vjDGbK5qYSJys4ikiUjakSNHziBWzZDQg5SEtmHPtk1WR1FKqRp10qJhjLEbY5oaY6KMMQ7/87LXzpPNKyILRWRDFY/xldZhqGKvRUS6AD2B9kA74DwRGXaCnC/5+8RKjouLq+YtB19kQhiIje/f/9DqKEopVaMCvU7jlBljxpxonIgcEpE2xpifRaQNcLiKya4AvjPG5Pvn+Qw4G1galMA1KPnycXz2dAZ5PxVYHUUppWpUoFeE17SP+eWCwanAR1VMswcYISIOEXHiOwheZfNUXdO5d19Cin/GFB/X6qaUUvWaVUVjJnC+iGwHxvhfIyLJIvKKf5r3gJ3AemAtsNYY84kVYU+HzezB7UzA7dIuRZRSDUfQmqdOxhiTCYyuYngacKP/uQe4pZaj1Rh7TB7FJZGs+OQDhl05yeo4SilVI6za02jw2iT5el75aflai5MopVTN0aIRJMMnXYvNU4wrU/ufUko1HJY0TzUG4U2aEFKyGyN6kZ9SquHQPY0gkpCfKQ1tx+G9u62OopRSNUKLRhCFt7NjbHaWzX/X6ihKKVUjtGgEUf+xF4Lxkr0p1+ooSilVI7RoBFGP5KGEFW3H4+mj12sopRoELRpB5mi+l9LQOBa88KzVUZRS6oxp0QiyEbdci3hdZPyoTVRKqfpPi0aQJfToQ2jJBly2vhTk5lgdRymlzogWjVoQFp+L29mUT598yuooSil1RrRo1IJLbr8Vu7uQgp2hVkdRSqkzokWjFsTGtcbpXktJSCIHd++0Oo5SSp02LRq1JKa3A689lC+ffc3qKEopddq0aNSSsbfdjrM0i9KDLa2OopRSp02LRi0JCQvDLuspCevJllXfWx1HKaVOixaNWtT+3PYYm53vZld1d1ullKr7tGjUogtSbyS0eD+e3E5WR1FKqdOiRaMWiQi2sC0Uh5/Fik8+sDqOUkqdMi0atazX+CEAbP4kzeIkSil16rRo1LKhF48jrGgHntJe2vOtUqre0aJhAXv0T5SGteF/s1+xOopSSp0SLRoWOCf1CsTr5uCKQ1ZHUUrVMo/X8MiCzdw8J40vNh7E5fFaHemUOKwO0Bh1GzCYpcX/xu3oS1FBAeFNmlgdSSlVC4pdHma8s4bPNhwkJsLJ/zYdokVkKBMGtGNicjxdWkZaHbFauqdhkdA2R3CFxPDp09rzrVKNQV6xi9TXVvLZhoP8+bJepN03hleuS6Z/hxheWfYTY574mqueX867aXspKHFbHfeExBhjdYYalZycbNLS6v6ZSYf37ub9v24kpHQtN7x+r9VxlFJBdCSvhNTXVrL1YB6PT+zH5f3bHTP+cF4x//lxP+/+sJddGQU0CbEztl9bJg2Kp398DCIS9IwissoYk1zddNo8ZZGW8R0JKX0bl7MvmQf307x1u+pnUkrVO3syC/n1rO85nFvCK1OTGdn9+P7nWkaFMX3EWdwyvDNpu4/yzg97+WjNAd7+YS+dWjSheZOQ8mkr/8yv+MO/W6soZk7oG6y3AmjzlKWiurnwOML5/KkXrI6ilDqBDftzWLz18GkdsN54IIcrn19OTpGLt24aUmXBqEhEGJTQjMcn9uOHP41h5pWJdGrRhFCnrfwR5rQR7rSXPyJCHDQJ9T3CnPbTfZsB0z0NC1162228ecdiivc3szqKUqqSwlI3//h8K6+vSMcYaBkVyqTkeK4eFE98s4hq51+xM5Ob56QRFebg7ZvPpkvLqFNaf2Sog2sGd+CawR0Cmr60qJCCjIOntI7ToUXDQk2aRuM06ygOHcKO9avpktjf6khKKWD5jgzu/s869mYVMfXsjpzTpQXv/rCX55bs4NklOxjWNY5rB8czumcrnPbjG2w+3/Azv397DR2aRTBn2mDaxoSf0vpLiwo5tG0dh3dsJ3v/IQqPFlKc58JVJLhLnHg9oXi94Rgi8Nqa4LE3wWMPJ6z4J254vXNNbYYqadGwWNygWPasc/LtK/Pp8m8tGkpZKa/YxSMLtjBv5R4Smkfw7i1nM7iTryXgwt6tOZBdxLtpe3n3h71Mf/NHWkSGMjG5PdcMiqdjc9+p8/NW7uG+D9bTLz6GWVMHEes/HuEuKeHg1rUc2LSRo3sPkZ9RSEkBuIudeNzhGNMEr0ThsUfidkSC2IB2/scvbFKMnQJsFGCjEKc5Sqgpxi4uIuJr4YC5nj1lLbfLxWs3v4fNk8MNc6ZbHUepRmvJ1sPc+5/1HMot5oZzO/GH87sTHlL1MQKP1/DNtiO8tXIPX205TLgrl/HRh2lbnEXB4TxiTAihJgKvOxKvicJji8HtbIqxHf873eYpxuHOw+bNw0Y+NnsRjpASnBEewmOcRDaPJLp1S5p16EjLLj2IiI0LyvvXs6fqCYfTid2xiaKQEaxa9BkDR19sdSSlGpWcQhcP/XcT763aR5eWkbx36zkM6BB7zDQZP20hPW0lGbv2kXe4mOI8O+7SSIaaGAbZY3E5W0Ce7yB3Exu4AK+7CAc52MghhF2ESxGh4W4imoUQ3bYZcZ0SaNMrieh6duakFo06oOvFiaz7n4117y3VoqFULfrfxoPc9+EGsvKKmNGtmH6unex+4Ss2Zwmu4kg8phkuRws8zkigvf8BggunHMVujhJittHEXkRYtMHRPIKEnl2I79ufZvFnWfregkWLRh0w7MpJbP1oFh7pRWlxMSFhYVZHUqpBSt+ygbVffkn2riO4chzgjuEWaYErJA7vyuZspLlvQuMlRLKwmwzCWU+Is4QmLRw06xBH+z59aN83BUdoqLVvxiJaNOqIkLid5OWP5u0/PcR1jz9sdRyl6i23y8WGb79m+4qV5O/Pw1sQAZ4WeBytcIXEAP0AEIcbpzcDu8mkidlDWHgJ0e2a0K5PN84aOpLwGD0VvipaNOqIa/72Z+bc8hbFpQP5OX0HbRK6WB1JqTpvz7ZN/Pjfz8jelYk3PwLjbY3L2QaPIxwYCoDdXoTDewi7ZxsO8ghr6aRd3+4kX3AJTZpGW/sG6iEtGnVESFgYLfpns39LAp89/DrTXn7I6khK1RmlxcV8t+AD9qZtouQImJLmeO1tKQ1tAfhOVbfbC3F69uN0ryY0rJDIdlF0GZpM4rkX4HA6rX0DDYieclvHzJr6OCUhfRg40c3gCy6zOo5Sta60uJjlH7/P3pWbcGWFYDxtKA2Jx2v3H0MwXkJKDmMzPyOhWYS3CeGslGQGjNLicCbq9Cm3IjIReADoCQw2xlT5LS8iFwH/BuzAK8aYmbUW0iJ9Jvdg1XxYP3ezFg3V4LldLlYu+Jhdy9dQesSG8bTB5YzH42gDtMFmLyXEs5dQ9w/YI4tp3r0tQ8aNJ67tGKujN1qW7GmISE/AC7wI3FlV0RARO7ANOB/YB/wA/MoYs+lky67vexoAs278M0WOUbTrsYbL7/iD1XGUqjF7tm3iu/kfkJ9ehClpjSukIx6H70pq8boIKdmPTQ5gjy2kTVJnUq6YqMcdakmd3tMwxmwGqusjfjCwwxizyz/t28B44KRFoyG46L7r+PSh9WSsjtFTcFW95Xa5WP7Jf0j/dh2uzAgM8ZSEtgU5G2xeQm0/E+Jeh71JAXF94km5ciLRzS+0OraqRl0+EN4O2Fvh9T5gSFUTisjNwM0AHToE1iNkXda2U1fCYmeTV6Cn4Kr6I/doNl/NeY2szRmYwpa4nQm4nXHAaOyOIpylPxHu3U5kQhhDJ15Bh27axFQfBa1oiMhCoHUVo+4zxnxUk+syxrwEvAS+5qmaXLZVrnn4l1NwD+7eSeuODfPqUlV/HT1ykCVvvMnRLdmY4raUhnbCa/ddAxFiO4TTs5mwyFzaJndhxMRrdY+5gQha0TDGnOnPiP1AfIXX7f3DGoWKp+Au+NtsPQVXWS7z4H4Wz3mD3O0FmNJ2lIYk4LUPABuE2vYT4k4jLM5F0rgL6T30V1bHVUFSl5unfgC6ikgnfMXiGuBaayPVrsvv+IP/FNwUVv7vUz2bqgHyeg1eY3BUcU8GqxUVFLBw9iscWXsYUxxPSVgnjG2w73gE+wn1fE9YKy/JV15KtwG/tjquqiVWnXJ7BfA0EAf8V0TWGGMuFJG2+E6tvcQY4xaR3wFf4DvldpYxZqMVea1UdgruBj0Ft8Hxeg03v5HG+v05PHPtAAYlWNtthdvlYvnH77FryXo8uXG4nF3xOBJ9exKyl1D3d0TEw6CrxtElUY9HNFZ6cV89oKfgNkyvLvuJhz7dREyEk/xiN/de0pNpKQnVnVVYozavXEba+59TerAJHltXXKG+wuUsycJuthPaqoCBEy6m5+CUWsukrFGnT7lVp0ZPwW14Nh7I4dHPtjCmZ0v+OSmJO+ev5aFPN/HjnqM8OqEvkaHB+a9ZkJvDl7NeIWtDNl7XWZSEdwDOw+4sxFm6nRDnj5w1Komzx16hV1erKumeRj0x5//dR17BaKIiF+kpuPVcUamHsc8sI7fIxed3DKdZkxC8XsOL3+zisS+20DkukhemDKBLy6gaWd+6ZYtZ8/6XuLKa43J2811MZ7yEFf+EhO4mrl8cY1JvJLxJkxpZn6qfAt3T0KJRT5QWFzPnlrfw2mMY95d+egpuPXbfB+uZ+/0e3rxhCOd2bXHMuOU7M/j9vNUUlXp49Kq+XNa37SkvvyA3hy9eeYnsjXl43WdREu47CdFZmo3ds5XQ1oUMnXIlXRL1nvTqF9o81cDoKbgNwxcbDzL3+z3cMrzzcQUD4JyzWvDpbcP4zdxV/O6t1fy4O5t7L+mBs5qzq3asTWP5mx/gPhxLqbMHHsdAcHgIc+0iwiyi9eAOnH/dNBzOK4P11lQjoXsa9UxZL7hJV5Ry9iXjrI6jTsHBnGIu+vc3tI8N5z+3phDiOHEhKHV7eWTBZmYvTye5YyzPTh5Aq6a/HMtyu1x88948dn+9BW9hR4rDO4PYcbhycbg3E9auiJTrJpHQo09tvDXVAGjzVAO18n+fsmq+Dacri3bXtuXiC86zOpIKgNdrmPLq96zek81/f38uneMiA5rv47UHuOf9dUSEOHjsko5kfD6fvO1uPHSnNDQOgNCivdicO2mR1JyLbrxFT5RQp0WbpxqowRdcxqEN/2Dvln7sn5vOP7MWcvvE8+rkxWHqFy8t3cXynZk8OiEx4IIBMCQ8kwfcy8lLb8rW9BI8jqHYHKWElGyjiXMNvccNY9D5U4OYXKljadGoh8b+4f/49O8Pszt9MDELdjDtcDGPpY45pvlC1R3r9mXz+BdbuSSxNZOS46udfsvi/7L24+/Iz2pHcdhZICNxOHIIca8lvF0p592YSqv4i2ohuVLH0+apeuzDPz/A/sPnEla0kVc7JfCPKSlVHlxV1ikocXPpU0spdXv57PbhREccf+2Dp9TFyndeY+e3ByguPouSsHYAhBb/TFjodhKGtmTotTfgCA2t7fiqEdHmqUbg8oceYP7/+zOHZRQ3pKcx9RUvvxvTm9+P7ordVntXFasTe+DjjezOKuTtm4YeUzCKsjNZ9uqrHNjkpsT0wBXSBUxnwsxOmocups/FA+hz0WQLkytVNS0a9dzEfz7EW7+7n6PhI/njoeU89KWdtN1ZPHl1f+Ki9JeplT5dd4D5q/Zx23ldGNK5OUd2bmb56++StSeSEkcPPI5kbLYSwkq3ENdsLUOmjKNtr1usjq3USWnzVAPgcbt56zcPk2sbRhPPEv7e6myahjt56pr+nH1Wc6vjNUr7jhZy8b+XMsrxEylH91OQ0Ybi0K4Ym++02DCzmVY9DOekptK05alfwKdUTdNTbhsZj9vNmzf/g/yQoUTbv2JW5zGkZxRw54XduXXEWbXaCV5j5na5WPTmbHYv+wmbuysl4R0BCCk+RFjINhIGN2fI5GmEhEdYnFSpY2nRaIRKiwqZe+tzFIYNIC5qMUsHXMknaw9w3dkdeWBsb2x6nCMo8rKz+fz558jb5sZDD0pDfScjhBX9RETUbnpe2JuksVdbnFKpk9MD4Y1QSHgEv3r6Rub97g2OmBGM2fIxbYZP4qVvdpFd6OLxif1OehWyCtyujev4ds67uA5GU+rsWeH6ia1Ehq6h/6Tz6Ztyg9UxlapxuqfRAOUePsD8uz6hJLQznTr/wO6zr2bmZ1sY0S2O56cMICJEfyucjuWffMDWL1bize1ASVgX//GJPBzuTYR3KGH0TVNp1aGT1TGVOi3aPNXIZaRv46O/fEtxWDyRrpU4xg3iz9+XkhQfw6zUQcREhFgdsc7LPXqU/730ArlbS/B4u1Ia1gaAkOKfsdu30axvNBfeeIt2Ka4aBC0aiqw9O/j0obfIs52Nw11IRMxK/hYylI4tmzJn2hBaR+sV5JWt+3Yxq9/7EndmC1wh3fE4whGvi9DiHdgi99NlTF+GXXmN1TGVqnFaNFS5H96dzbr/uikO70xY0U5WNfuZda3P5Y0bhtCpReP+lVxUUMCXs14mY10m3tKyO9n57j3h8G4mrH0pI2+YQttOXS1OqlRwadFQx3CXlPDxXx7myJFk3I5wnCXLeD++G/+cfhl92kVbHa/WuF0uvv1oPj99vQFPbhwuZ1c8jogKd7JLp82QsntP6O1OVeOhRUNV6efNq/nfY5+THzIEZ+lRskJXMHzGXZzTrZXV0YJm3beLWfPhQlyHI/HYu+IKaQaAsyQTu9lBSKsChlw9jm4DBlucVCnraNFQJ7X01WfZujSKkrD2hBVuotn5MYyfPJmMghIO5ZRwMLeYgzlF/n9LOJRbzMHcYg7lFBMd4eSW4Z2ZmBxPmNNu9Vs5zrrlX7NhwVcUH7Dh9XYu7wDQ7srH6d6OPTqDLqMHMPSSy3VvQik/LRqqWiX5efznT/8gJ38oXpuD0OKNHA3dzxdNz2K703cls90mxEWG0io6jNZNQ2ndNIwNB3JZtfsorZqGMn3EWfxqcAfLikdedjbfvD2Xwxv24c1rhtfWsfziOpunlJCSXUjEflr1b8t5U1L1TCelTkCLhgrY5m+/YsUry3FLL1whMf72/XTCInfT/bzu9B9/DXbHL9d2GGNYsTOTfy/azvc/ZdEiMpTpIzpz7ZAOQb0GxO1ysen7ZWz4YgnFB8C42lIa2hGv3Xf6sLM0G7s7HWmSQbPuLRh2zdU0b90uaHmUaki0aKhT5i4pIW3+HHZ+u4+iwk7lZxKFlGQQat9Gm96hnJt6PeGxv9yz47tdmTz91Xa+3ZFJ8yYh3DS8M78e2pEmoadfPApyc/hx4WfsW7uF4kOleIuigBa4HK3wOH13vROvm5CSfdhs+whpWUq30UMZMOpCbW5S6jRp0VBnbNd3i/nx/cXkHm5BSUg3vPYQbJ5iQkvTsdnysTuKcYS5CIuy4Q53sq1ESMsP53BUR649rx/Xnd2RqDDfl3hRQQEH0rdzZPducg4doiAzh5K8QtwFLjzFXrzFDiiNxWtrSWlIHMb2S9FxlObgcB9C7BlIeCGxXZtz7jVX06JNe6s2jVINjhYNVaNyft7HijlzOLRVcHta4bFF4XFElTcNVWZ3F2L3FGDEgccejtdezYWExkNIyRFs3sOI8yj2pi5iOseReN5oOvfuG4R3pJSqSDssVDUquk17Lrr7j8cM87jd5P68h5+3biBz917yDmdTdLSEkgIodTnxesMwuHFQCrgRpwdbKDjC7YREhhMWE0XTls2JbduOTr0GEBUTY82bU0oFTIuGOm12h4PY+M7Exne2OopSqpZoP9lKKaUCpkVDKaVUwLRoKKWUCpgWDaWUUgHToqGUUipgWjSUUkoFTIuGUkqpgGnRUEopFbAG142IiBwBdp/BIloAGTUUJ9jqU1aoX3nrU1aoX3nrU1aoX3nPJGtHY0xcdRM1uKJxpkQkLZD+V+qC+pQV6lfe+pQV6lfe+pQV6lfe2siqzVNKKaUCpkVDKaVUwLRoHO8lqwOcgvqUFepX3vqUFepX3vqUFepX3qBn1WMaSimlAqZ7GkoppQKmRUMppVTAGk3REJGLRGSriOwQkXuqGP8HEdkkIutEZJGIdKwwziMia/yPj+tI3lQROVIh140Vxk0Vke3+x9Q6kPVfFXJuE5HsCuNqdduKyCwROSwiG04wXkTkKf97WSciAyqMq9XtGmDeyf6c60VkuYj0qzAu3T98jYgE/R7IAWQdKSI5Ff7e91cYd9LPkAVZ76qQc4P/c9rMP65Wt6t/nfEistj/HbVRRG6vYpra+ewaYxr8A7ADO4HOQAiwFuhVaZpRQIT/+a3AOxXG5dfBvKnAM1XM2wzY5f831v881sqslaa/DZhl4bYdDgwANpxg/CXAZ4AAQ4Hvrdiup5D3nLIcwMVlef2v04EWdWjbjgQ+PdPPUG1krTTtWOArq7arf51tgAH+51HAtiq+E2rls9tY9jQGAzuMMbuMMaXA28D4ihMYYxYbYwr9L78D2tdyxoqqzXsSFwJfGmOyjDFHgS+Bi4KUE04966+AeUHMc1LGmG+ArJNMMh6YY3y+A2JEpA21v10DymuMWe7PAxZ/bgPYtidyJp/303KKWS39zAIYY342xvzof54HbAbaVZqsVj67jaVotAP2Vni9j+M3eEU34KvYZcJEJE1EvhORy4MRsJJA807w74a+JyLxpzhvTQl4ff4mv07AVxUG1/a2rc6J3k9tb9fTUflza4D/icgqEbnZokyVnS0ia0XkMxHp7R9WZ7etiETg+4J9v8JgS7eriCQA/YHvK42qlc+u43RnbKhEZAqQDIyoMLijMWa/iHQGvhKR9caYndYkLPcJMM8YUyIitwCvA+dZnKk61wDvGWM8FYbVxW1b74jIKHxF49wKg8/1b9uWwJcissX/C9sqP+L7e+eLyCXAh0BXC/MEYizwrTGm4l6JZdtVRCLxFbA7jDG5tbHOyhrLnsZ+IL7C6/b+YccQkTHAfcA4Y0xJ2XBjzH7/v7uAJfiqfDBVm9cYk1kh4yvAwEDnrWGnsr5rqLSbb8G2rc6J3k9tb9eAiUhffJ+B8caYzLLhFbbtYeADfM1AljHG5Bpj8v3PFwBOEWlBHd62nPwzW6vbVUSc+ArGXGPMf6qYpHY+u7V5MMeqB749ql34mkbKDrT1rjRNf3wH47pWGh4LhPqftwC2E/yDdIHkbVPh+RXAd+aXg14/+XPH+p83szKrf7oe+A4gipXb1r+uBE58sPZSjj2YuNKK7XoKeTsAO4BzKg1vAkRVeL4cuMjirK3L/v74vmj3+LdzQJ+h2szqHx+N77hHkzqwXQWYAzx5kmlq5bPbKJqnjDFuEfkd8AW+MzVmGWM2ishfgTRjzMfAY0AkMF9EAPYYY8YBPYEXRcSLb89spjFmUx3I+3sRGQe48X2wU/3zZonIQ8AP/sX91Ry7a21FVvD9Ynvb+D/FfrW+bUVkHr6zeFqIyD7gL4DT/15eABbgOwtlB1AIXO8fV6vb9RTy3g80B57zf27dxtfLaSvgA/8wB/CWMeZzi7NeBdwqIm6gCLjG/3mo8jNkcVbw/Rj7nzGmoMKstb5d/VKAXwPrRWSNf9gf8f1oqNXPrnYjopRSKmCN5ZiGUkqpGqBFQymlVMC0aCillAqYFg2llFIB06KhlFIqYFo0VJ0nIvkBTHOHv8uHmlrn5SLSqwaXt/wM5s33/9tWRN47yXQxIvKb012PUoHQoqEaijuAUyoaImI/yejLgRorGsaYc2pgGQeMMVedZJIYQIuGCiotGqre8N+PYYm/g8YtIjLXfw+B3wNtgcUistg/7QUiskJEfhSR+f4+e8ruhfCoiPwITBSRm0TkB38neu+LSISInAOMAx7z3zPhLBFJ8nequE5EPhCRWP/ylojvfiFpIrJZRAaJyH/89y34W4Xs+RWe3y2++zGsFZGZVbzPTv7s6ystI0H8938Qkf/f3tmFWFVFcfz398FpaExhCIMRHPADwUGEehCkkUCEniyKIB9y1NUAAANKSURBVCY1fRWZl6wHKUjElHnTp/RBJdS+KBoi1IIwGoYZSvHOiEEx9SAJisjQSyLO6mGtq2eO98qZSRmcWT/YsDh773XWOffCumfve/5rtaThiK8maQVwAFgWx/oktclrw1wIX5sLfq5IOiqvzXBOUmv0LZf0Q8R2QdKyOL477lNN0oeP9INNniwe9+vv2bL930bU3MDf4B3HtXPmAYO4eBwUahzgkiQ/EfIPwHvAB4Vx7xZ8txfsfcCusI8Drxf6asCGsPcScg64XtbBsHuBv/HaBy24mmh76RpexqUn6rVbHpBzAPqBrWHvLMztJGQvgMNAT9jzgVZKshj4G8vPFO7JH7jERCeuJLA2+j4H3gp7CHg17Kfwp7dNwJGYOw/4Fuie6e9Ftplpc0JGJJlVDJvZVYCQU+gEfi6NWYcvLQ2E3MN8PMHU+axgd8Wv+UW4jMzZ8gklLQQWmdn5OHQC+KIwpC6VMgJcNrNrMW8MF4q7WRi7EThmUbvFGss5rAdeC/sT4GCDMYPAHklLgK/M7Pe41kmhA/sldQMTuBz24uj708zqchS/Ap2SFgAdZvZ1xPZvXMcmPHFcjPFtuDrtTCrmJjNEJo3kSeN2wb5L4++w8KIzbzbxUdQSOg68YmaXJL2NP81MN6aJUnwTTeKrwkP1fczslKQhXKTuO7k8/lhpWA/wLPC8md2R9Bf+9FCMGfw+tj7kdAI+MrOPpxB/MkvJPY1ktvAPXgYTvILdeknLASQ9LWllk3kLgGty2emeRv7MbBy4JenF6NsCnGd6fA9sr//TS1F3usQALvBIKaZ7yOuPjJnZIeAbYA2T7wG4Suv1SBgvAUsf9HQf84pwVxXFsCS1RJxngR2FfaEOeS2JZA6SSSOZLRwBzkj60cxu4Kq/pyXV8KWcVU3mvY+v4w8AvxWOfwrslnQxNoO34RvjNWAtvq8xZcwVUfuBX2J57Z0Gw3qBnZJGaF5h7Q1gNHx04WU+b+JLcqOS+oCTwAvhZ2vp+pqxBVdQruF7L8+Z2TngFDAYvr5kcnJK5hCpcpskSZJUJp80kiRJkspk0kiSJEkqk0kjSZIkqUwmjSRJkqQymTSSJEmSymTSSJIkSSqTSSNJkiSpzH/jSrr/3FkMJQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fed5ff03160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "for value, bootstrap in conditions.items():\n",
    "    plt.plot(results[f\"points_{bootstrap}\"], results[f\"energies_{bootstrap}\"], label=f\"{bootstrap}\")\n",
    "plt.plot(results[\"points_np\"], results[\"energies_np\"], label=\"numpy\")\n",
    "for condition in extrapolators.keys():\n",
    "    print(condition)\n",
    "    plt.plot(results[f\"points_{condition}\"], results[f\"energies_{condition}\"], label=condition)\n",
    "plt.legend()\n",
    "plt.title(\"Dissociation profile\")\n",
    "plt.xlabel(\"Interatomic distance\")\n",
    "plt.ylabel(\"Energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare number of evaluations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no bootstrapping\n",
      "Total evaluations for no bootstrapping:\n",
      "21482\n",
      "bootstrapping\n",
      "Total evaluations for bootstrapping:\n",
      "3908\n",
      "np\n",
      "Total evaluations for np:\n",
      "0\n",
      "poly\n",
      "Total evaluations for poly:\n",
      "1411\n",
      "diff_model\n",
      "Total evaluations for diff_model:\n",
      "3351\n"
     ]
    }
   ],
   "source": [
    "for condition, result_full in results_full.items():\n",
    "    print(condition)\n",
    "    print(\"Total evaluations for \" + condition + \":\")\n",
    "    sum = 0\n",
    "    for key in results_full[condition].keys():\n",
    "        if condition != \"np\":\n",
    "            sum += result_full[key].raw_result.cost_function_evals\n",
    "        else:\n",
    "            sum = 0\n",
    "    print(sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td><code>qiskit-terra</code></td><td>0.19.0.dev0+cdea7de</td></tr><tr><td><code>qiskit-aer</code></td><td>0.8.1</td></tr><tr><td><code>qiskit-ignis</code></td><td>0.5.1</td></tr><tr><td><code>qiskit-aqua</code></td><td>0.10.0.dev0+649b495</td></tr><tr><td><code>qiskit-nature</code></td><td>0.2.0</td></tr><tr><td><code>qiskit-finance</code></td><td>0.3.0</td></tr><tr><td><code>qiskit-optimization</code></td><td>0.3.0</td></tr><tr><td><code>qiskit-machine-learning</code></td><td>0.3.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.6.2 |Continuum Analytics, Inc.| (default, Jul 20 2017, 13:51:32) \n",
       "[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>1</td></tr><tr><td>Memory (Gb)</td><td>5.827335357666016</td></tr><tr><td colspan='2'>Fri Aug 13 14:31:56 2021 EDT</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2021.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
